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Sloshing of Liquids in Circular Canals and 
Spherical Tanks' 


BERNARD BUDIANSK Y* 
Harvard Unwersity 


Summary 


Theoretical calculations are made of the natural modes and 
frequencies of small-amplitude sloshing of liquids in partially filled 
circular canals and spherical tanks. An integral-equation ap- 
proach is used to analyze the circular canal for arbitrary depth of 
liquid. A similar approach for the spherical tank provides modes 
and frequencies for the nearly full and half full cases. These 
results, together with the known behavior of the nearly empty 
tank, are used in conjunction with the trends established for the 
circular canal as a basis for estimating frequencies for arbitrary 
depth of liquid in the spherical tank. The dynamic analysis of 
the container-fluid system by means of the mode-superposition 
approach is discussed, and modal parameters required in such 
analyses are presented 


Symbols 


geometrical coordinates (see Fig. 1 

nondimensional coordinates 

canal or sphere radius 

geometrical parameters related to liquid depth 

velocity distributions in x, y, and z directions, 
respectively 

oscillation frequency 

time 

velocity potentials 

displacement potential 

acceleration of free fall relative to container 
wall 

frequency parameter: \ = w?R/g 

liquid depth (variable) 

complex variables 

complex potential 

kernel function (circular canal) 

kernel function (spherical tank) 

mode shape for circular canal 

mode shape for spherical tank 


Received February 2, 1959. 

} This paper was prepared under Air Force Contract No. 
AF-04(647)-181 during the summer of 1958 while the author was 
a consultant to the Missile Systems Division of the Lockheed 
Aircraft Corporation. 

* Associate Professor of Structural Mechanics. 


= WV ‘pf p) 
= cylindrical coordinate system 
nondimensional coordinates 
p/p 
pp 
complete elliptic integrals of first and second 
kind, respectively 
Neumann function for sphere 
an, Bn modal parameters 
AB.C.D. 


M. mass of container 


nondimensional modal parameters 


PL density of liquid 

M, mass of liquid 

X external transverse force acting on container 
Fs sloshing force acting on container 

U transverse displacement of container 

A square matrix 


ig integrating matrix 
i column matrix 
| 


| row matrix 
integers 
A = interval of subdivision (A = 1/N) 


(1) Introduction 


KF R sufficiently small amplitudes of motion, the dy- 
namic effects of the sloshing of a liquid in a par- 
tially filled container can be analyzed in terms of the 
natural modes and frequencies of the small, free-surface 
oscillations of the liquid. This paper is concerned with 
the calculation of the modes and frequencies of par- 
tially filled circular canals and spherical tanks, with 
attention restricted to the modes that would be induced 
by lateral accelerations of the containers. The usual 
assumptions of inviscid, incompressible flow are made, 
and the container is considered to be rigid. The dy- 
namic analysis of the container-fluid system by means 
of the mode-superposition approach is discussed, and 
modal parameters required in such analyses are pre- 
sented. 

The case of the half full circular canal was considered 
by Rayleigh (see reference 1, p. 444), who carried out 





162 JOURNAL OF 





FREE SURFACE F 





WETTED SURFACE S 





Fic. 1. Container partially filled with liquid. 


an approximate energy solution for the fundamental 
frequency. Graham? estimated the fundamental fre- 
quency of the half full sphere by conceptually dividing 
the liquid into a collection of half-discs and then using 
Rayleigh’s result in an unspecified averaging process. 

The frequency and mode calculations of this paper 
are based on a rigorous integral equation approach. 
In the case of the circular canal, the appropriate in- 
tegral equation is formulated explicitly for an arbitrary 
depth of liquid, and numerical calculations of the first 
three pertinent modes and frequencies are performed 
for a range of depths by means of a matrix approxima- 
tion to the integral equation. For the spherical tank, 
rigorous analyses are made for only the nearly empty, 
half full, and nearly full tanks; results for the inter- 
mediate depths are estimated by interpolation and by 
means of the trends established for the circular canal. 

The results for the fundamental frequency for the 
spherical tank are compared with available experimental 
data. 


(2) Sloshing Modes and Frequencies 


(2.1) Basic Differential Equations 
The differential equation and boundary conditions 
governing free-surface oscillations (reference 1, p. 363) 
will be displayed. 
The liquid is supposed to be bounded by a horizontal 
free surface F and a rigid wetted surface S (see Fig. 1), 
and small, harmonic oscillations are assumed. The 
velocities in the x, y, 2 directions, respectively, are then 
v,(x, ¥, 2) sin wt 
v,(x, ¥, 2) sin wt (1) 
v(x, ¥, 2) sin wt 

where w is the circular frequency. 

Since the flow is irrotational, there exists a velocity 

potential 
o(x, y, Z) sin wt 
such that 
Uv, = 0¢/Ox, 


v, = 06/dy, v, = 06/02 (2) 


The continuity equation for incompressible flow is 
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V4 = 0 (3) 


The condition of tangential flow at the rigid surface 
requires that 


Ag-n = 0¢/dn = Oon S (4) 


where n is the unit normal vector to S. 

Finally, the condition of constant pressure at the 
free surface together with the linearized Bernoulli 
equation leads to the condition 

0¢/Ov = (w?/g)d on F (5) 
where g is the acceleration of free fall relative to the 
container. 

Eqs. (3), (4), and (5) govern the eigenvalue problem 
for the frequency w. Note that Eqs. (1), (2), and (5) 
imply that the shape of the vertically displaced free 
surface during oscillation is given by ¢(x, 0, 2). 

The analogous two-dimensional problem of trans- 
verse sloshing in a uniform channel is formulated in an 
obvious fashion by suppressing the z coordinate. 


(2.2) Circular Canal 

The geometrical parameters for the partially filled 
circular canal are shown in Fig. 2; the depth of liquid 
is measured by e, which varies from —1 for the empty 
canal to 1 for the full vessel. 

It is clear that, because of geometrical symmetry, 
the natural modes involve free-surface vertical dis- 
placements that are either symmetrical or antisymmet- 
rical about the origin. Since the symmetrical modes 
would not be induced by transverse motion of the canal, 
attention is restricted to the antisymmetrical modes. 

It is convenient to nondimensionalize the problem by 


introducing 
& = x/aR, » = y/aR (6) 


The normalized geometry is shown in Fig. 3; the 
free-surface condition becomes 











0¢/O0n = (Aa)d (7) 
where \ = w?R/g. The angle 6 shown in Fig. 3 is given 
by 

y 
oR 
_——$$_$_—__—_— ar} 
\ on 
R 
Fic. 2. Geometrical parameters and coordinates for circular 


canal and spherical tank. 
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Fic. 3. Nondimensional geometrical parameters and coordinates 
for circular canal and spherical tank. 


B = arc sin e (8) 


and varies from 7/2 for the full canal to — 7/2 for the 
empty canal. 

The procedure for setting up an integral equation 
governing the antisymmetrical modes follows. 

Consider the fluid-filled container consisting of the 
original wetted boundary plus its reflection about the 
t axis (Fig. 4). Place a two-dimensional sink of 
strength 2w(=)d= at E and a source of equal strength at 
—£ and smear out both over a strip of length d=. De- 
note the resulting velocity potential for the internal 
flow by 


w(E)@(E, ; E)dE 

Since the velocity potential is arbitrary to within a 
constant, it will be assumed that ®(0, 7; —) = 0 and, 
hence, ®(£, 7; 2) = —(—E, n; &). By symmetry, the 
vertical velocities along » = O vanish everywhere ex- 
cept on the strips of length dé at — and —é, where just 
below the £ axis they are w(E) and —w(é), respectively. 
It is now possible to construct the velocity potential 
for the free-surface oscillation of the configuration of 
Fig. 3 by means of a distribution of sources and sinks 
along the £ axis. It is necessary only to satisfy the free- 
surface condition, Eq. (7), relating the vertical velocities 
along the ~ axis to the velocity potential along the & 
axis. Thus, we require that 


1 
w(t) = (Aa) f &(é, 0; E)w(E\dE (9) 

0 
for0 < — < 1; since (—£, 0; —) = —(E, 0; &), Eq. (9) 
is automatically satisfied for negative £, with w(—£) = 


—w(§). 


The problem now is to find the kernel function 


G(é, £) = #(E, 0; &) (10) 


and this is readily accomplished by conformal mapping, 
The successive transformations shown in Fig. 5 map 
the region of Fig. 4 into the entire f plane; the re- 
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sultant transformation is 
f = {(01/2)(Q — 2)/ + 2)]}? 


where y = 2/[(1 + 28)/z]. 
spond to z = & and —&, respectively. 
potential of a sink at f; and a source at fo, each of 


(11) 


Now let f; and f2 corre- 
The complex 


strength 2, is 
F = —(1/m) log [(f — fi)/(f — fe)] + g(&) (12) 
where g is an arbitrary function of £. 
Substituting Eq. (11) into Eq. (1 
plane, the corresponding potential 


2) gives, in the z 


] 
F=- log 
TT 
(a — 2/1 +217 - (1 -5)/0+5l + of 
} ef 5 - en % y(E) 
ia —2/a+s2))7 —{a+8)/a -— Bs & 


Now choose g so that F(0) = 0; then 


» t 


Fic. 4. Location of source-sink pair in symmetrized domain. 
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Fic. 5. 





164 JOURNAL OF THE AERO/ 


poneae 


a c 


. 6. Limiting case of nearly full container. 
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log 
rs 


y((1 — 2) 
({a — 2) 


+o -la- 
(1+ s))* — (1+ 
fl-— fa 
i -—[a—- 


Serl | Seri 
4. Lees 


| 
— log 
T 


zy"(1+é ied. 


ja - 
(1 + 2)"(1 ‘—s 


Re F(z) is 


z)"(1 — a 
z)7(1 — &)” 


Finally, G(é, —) = 


a+epai-—-)’-a-p7a+ 8’ 

aga+eaa+5?-(d-—é (1 — 

Note that G has a logarithmic singularity at & = &. 
Note also that 


(13) 


G(é, £) = G(é, §) = —G(—§, &) 
To recapitulate, then, the solutions of the homo- 
geneous (singular) integral equation 
1 
w(t) = (na) f G(é, E)w(E\dé (14) 
0 


with G(é, —) given by Eq. (13), provide the character- 
istic free-surface shapes w(¢), and the associated eigen- 
values of (Aa) yield the natural frequencies. 

The numerical solution of the integral equation, 
Eq. (14), has been performed for various depths on the 
basis of the matrix formulation described in the 
Appendix. The results are presented in Section (4). 
Special Cases—Several special cases for which approxi- 
mate solutions can easily be obtained are of particular 
interest. For the half full canal (e = 0), the kernel 
becomes 


|}€-—H0 - & 
= Jog | = 
(— + &)(1 + &€) | 
A one-term Galerkin-type solution for the lowest fre- 
quency can be carried out without difficulty for this 


case. 
From Eq. (14), 


1 1 
f [w(é)]? dé = (da) f f "G(& B w(® w(t) aedé 


(15) 


SPACE 
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Introducing the approximate mode shape w = ¢ 
and carrying out the integrations give, 
the fundamental 


into Eq. (15) 
ior ¢ = ©. fj = 1), 
eigenvalue 


approximate 


hi = 1/3[(2/r7) — (2/8)] = 1.367 (16 
or 
(g/R)'/? 


w, ~ 1.169 


This is precisely the result found by Rayleigh in an 





entirely different way. The result is necessarily an 
upper bound to the exact frequency. Rayleigh also 
obtained an improved approximation w ~ 1.1644 | 
(g/R)'”?. 

Also of special interest is the limiting case of the 
nearly full canal, where 8 ~ 7/2 and y > 1. The 
integral equation becomes simply 

1 | ‘eee: 
w(t) = — 0) f log — = w(t) d= (17) 
0 7 E+ ¢ 
Denote the critical values of (Aa) by c, (z = 1, 2,...); 


then, for e sufficiently close to 1, 
An ~ Cn/a 
= Ce V1 —@¢ 


Thus, the frequencies tend to infinity as the full 
This result, as well as the 
(17), could have been 


condition is approached. 
form of the integral equation, Eq. 
anticipated by consideration of the sloshing of liquid 
in a half plane with a rigid boundary having an aper- 
ture of length (2aR). (See Fig. 6.) The kernel func- 
tion for this case can be written directly in terms of a 
source-sink pair, without the need for conformal trans- 
formations. A one-term Galerkin solution of Eq. (17), 
again with w = &, gives the approximate result 

~ 2/3 = 


2.094 


so that 


2.094/VW1 — e (18) 


A =~ 


fore 1. 

For the nearly empty canal 
infinite. The asymptotic results for this case are most 
readily found by abandoning the integral equation and 
exploiting shallow-water theory (see reference 1, p. 


(e —~ —1), y becomes 


291). The differential equation for the velocity po- 

tential may be conveniently deduced from the vari- 

ational principle given by Lawrence, Wang, and 

Reddy.* The variational equation is 

i) lf (1/2)|V¢) 2dV — (w?/2g) f ¢°* aF| = (19) 
Vv F 


where V denotes the volume of fluid in Fig. 1. 

As shown in reference 3, Eqs. (3), (4), and (5) consti- 
tute the Euler equation and natural boundary condi- 
tions of Eq. (19). 

For the two-dimensional canal of shallow depth vari- 
ation h(x), it is assumed that @ may be taken inde- 
Then Eq. (19) becomes 


pendent of y. 
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if » f h(x)(do/dx)? dx — (w?/2g) i} dé? ax | = 0 
F I 


The Euler equation is 
(d/dx) |h(d¢/dx)| + (w?/g)=¢ 0 (20) 
For the nearly empty circular canal, the depth has 
approximately the parabolic variation 
h(x) = (1/2R)[(aR)? — x?| 
so that Eq. (20) becomes 
(1/2)(d/dx)[(aR)? — x?(do/dx)| + (R/g)o = 0 
or, in terms of § = x/aR, 
(1/2) (d/dé) [(& — 1)(db/dg)] — Ad = O 
whence 


(€? — 1)¢” + 28’ — 2 = 0 G1) 


This is essentially Legendre’s differential equation. 
For the antisymmetrical modes, (0) = 0; in addition, 
the condition of finiteness at & = 1 is imposed. The 
solutions of Eq. (21) are then the odd-order Legendre 
polynomials of the first kind 


Poi() =» = 1,2,3,... (22) 
with the corresponding eigenvalues 
L, = 2? — 2 w= 1,23... (23) 


The polynomials [Eq. (22)] represent the limiting 
free-surface shapes of the natural modes as the depth 
approaches zero. 

The first eigenvalue is \; = 1; and so the fundamental 


frequency 
Qy = V2/R 


is, as could be expected, simply the frequency of a 
pendulum of length R. 


(2.3) Spherical Tank 
The geometrical parameters of Figs. 2 and 3 may be 
considered to apply to the spherical tank. A rigorous 
analysis was executed for only the nearly full and half 
full cases; rigorous results for the nearly empty case, 
given by Lamb (reference 1, p. 291), are also available. 
Introduce the cylindrical coordinate system (r, y, 8), 
where x = r cos 6,2 = ry sin 6. The appropriate non- 
dimensional system becomes (p, 7, 6), where p = r/aR, 
and, as before, 7 = y/aR. In these nondimensional 
coordinates, the governing differential equation is: 
(0°6/Op?) + (1/p)(0b6/Op) + 
(0°p/On*) + (1/p*)(0*¢/06?) = 0 


and the free-surface boundary condition remains that 
given in Eq. (7). It is evident that solutions for the 
modes may be written in the form 


(k = 1,2,...) 


(24) 


(25) 


It will be shown later that only the modes with k = 1 
are induced by lateral accelerations in the & direction. 


& = ¥(p, n) cos k6, 
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i 








Fic. 7. Coordinates in free surface of spherical tank. 


Consequently, attention will be restricted to these 
cases. 

The procedure for setting up an integral equation 
analogous to Eq. (14) is as follows. Consider the vol- 
ume of fluid enclosed by the w~’ 
reflection about the originally free surface » = 0. 
(See Fig. 4.) Now, in the plane 7 = 0 (Fig. 7), place 
a distribution of three-dimensional sinks along the 
annulus (f, p + dp) of strengths 2/(f) cos @ per unit 
Let /(A)2(p, , 6; p)pdp be the resulting velocity 
potential. Then, since the vertical velocity just below 
the plane » = O is, because of this potential, zero 
everywhere except in the annulus, where it is /(j) cos 6, 


ed surface and its 


area. 


the imposition of the boundary condition Eq. (7) at 
6 = Orequires that 


1 
f(p) = (Aa) i} Q(p, O, O; p)pf(p) dp (26) 

0 

It is convenient to let g(p) = Vv» f(p); then, 
1 
g(p) = na) f H(p, B)g(b)dp (27) 
where 

H(p, B) = V ppQX(p, 0, 0; 3) (28) 


The problem now is to determine //(p, A). 

2.3.1) Nearly Full Tank—For the nearly full case, 
the symmetrized region is the entire space (Fig. 6), and 
the potential of a sink in the plane 7» = 0 is unmodified 
by the presence of a rigid boundary in this plane. The 
potential of a sink of strength 2 at (A, 0, a) is 

x 
2r 
] 
[(p cos 8 — pcos a)? + (psin 6 — psin a)? + 7?]'”” 
and so 
l cos a da 
Q(p, 0, 0; p) = P -f 9 ~s od 1/2 
2rJ) [p? + B? — 2pf cos a] 
Then, 
5 1 " cos ada 
Hp, p) = an V op -¢ ae 
T o lp 


24+ 5? — 2pp cos al'”” 
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Geometrical parameters in equatorial plane of sphere. 


Clearly, H(p, 6) = H(j, p). Consider p < §, and 











let m = (p/p). Then, 
ta, Vim . — cosada— - 
« Jo [m? — 2m cosa + 1)’ (29) 
Letting 8 = a/2 yields 
t/2 (9 on2 a> i 
H=41 f Sach wd (30) 
xJo [1 — g? sin? p]°’~ 


where g = 2\/m/(m + 1). 
This may be evaluated in terms of complete elliptic 
integrals (reference 4, p. 73) to give 


H = (1/m){[(2/¢) — J K(@ — (2/QE (g)} 


where 


(31) 


K ggueen...ae 
(q) = J [1 — g® sin? 6]'/* 


and 


19 


E(q) = f " [1 — q? sin? 6]? do 


are the complete elliptic integrals of the first and 
second kind, respectively. 

The result [Eq. (31)] for H(p, 6) may be put into a 
much more convenient form by exploiting the Landen 
transformations? 


K(m) = [1/(1 + m)]K(g) 
E(m) = [(1 + m)/2|E(q) + [(1 — m)/2]K(q) 
to give, asa final result 
H(p, B) = 2/V 0/BIK (0/2) — E(0/2)] 


for p <p. Since H(p, j) = H(j, p), the kernel is com- 
pletely defined. This kernel, like the ones for the 
circular canal, also has a logarithmic singularity at 
p = p. 

(2.3.2) Half Full Tank—For the half full case, the 
symmetrized volume is a complete sphere, and so the 
kernel function can be evaluated on the basis of the 
known Neumann function (or Dirichlet function of the 
second kind) for the sphere. Thus,® the function of 


(32) 


the position Q (Fig. 8) given by 
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N a fp oe + 
a = 2¢ \PQ °° (OP)(QP’) 
9 
~ ( 


 ——_ Fae oe Ot 
°8 1 — (0Q)(OP) cos y + (OP)(QP’)S \°?) 
satisfies the following conditions in the unit sphere: 


(2) N(Q) is harmonic. 
(771) N(Q) hasa sink of strength 2 at P. 





(1172) The normal derivative of V(Q) has the constant 


value —(1/27) at the surface of the sphere. 


Now, let Q be the point (p, 0, 0) in the equatorial 
plane of the sphere (Fig. 7), and let P be the point 
(p,0, a). Then, 


1 j 1 
a a Les 
. Qe U[p? + p? — 2pp cos a]/? 
] 


| 


Cee or 


2 


log 


1/2 


1 — pp cos a + [(pp)? + 1 — 2pf cos al! 
and, therefore, 


H(p, B) = 
V pb { cos ada 


» [p? + p? — 2p cos al '/? 
T / —_ rT 
cos ada V pp 
f a : a : cos a X 
0 [(pp)? + 1 — 2pp cos al!” T 0 


) 
tog ——— : | da 
1 — ppcos a + [(pp)? + 1 — 2p) cos al?/? 
(34) 


+ Ve 
Tv 


us 


Note that, because of the cosine variation of local 
sink strength at the radius j, the normal flow at the 


spherical surface vanishes, as it should, despite the 


fact that the Neumann function [Eg. (33)] for a single 
sink produces a constant flux at the surface. 





The evaluation of the first two integrals in Eq. (34) | 


follows immediately from the results for the nearly full 
case; the last integral may be evaluated as follows. 
Let s = (pp), and let 


I(s) = i) cos a X 
0 


log J — - — —- 
ue — scosa + (s? + 1 — 2s cos a)|?/? 


7 
--f cos a X 
0 


log [(1 — s cos a) + (s?+ 1 — 2s cos a)!/*Jda (35) 


9 


t da 


Note that 7(O) = 0; subtract the zero quantity 


f (cos a) log sda from Eg. (35), giving 
0 


I(s) = -f cos a log }(1/s) — cos a + 
0 
[1 + (1/s?) — (2 cos a/s)]/?} da 


Then, by differentiation under the integral sign, it is 
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f = Oat p = O, and f is finite at p = 1. 
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found that 
, ] { cos ada 
I"(s) = : 
s Jo [s? + 1 — 2s cos a]}/? 
which [see Eqs. (29) and (32)] is 
I(s) = (2/s*)[K(s) — E(s)] 


Hence, 
i(s) = { (2/s*)[K(s) — E(s)|ds 
0 


which gives (reference 7, p. 273) 
I(s) = (2/s)[E(s) — (1 — s*)K(s)] 
Then, with m = (p/p)(p < j), evaluation of the three 
integrals in Eq. (34) gives: 
H(p, B) = (1/m)} (2 Vm) x 
[K(m) — E(m)] + (2/Vs)[K(s) — E(s)] + 
(2/V's)[E(s) — (1 — s2)K(s)]} 
or, finally, 
™) Xx 
{(2/V p/p) [K(o/B) — E(0/p)] + 
2(pp)*/* K(pp)} 


H(p, 6) = (1 
(36) 


for p <p. Again, the additional fact that H(p, p) = 
H(p, p) [see Eq. (34)] completely defines /7(p, 5). 

As in the case of the circular canal, the solutions of 
the integral equation with the kernels Eqs. (32) and 
(36), were performed on the basis of the matrix setup 
discussed in the Appendix. The results are given in 
Section (4). 

(2.3.3) Nearly Empty Tank 
tank may be analyzed in a fashion analogous to the 
treatment of the nearly empty circular canal by ap- 
proximating the wetted surface with the paraboloid of 
revolution having the depth variation 


The nearly empty 


h(r) = (1/2R)[(aR)? — r?] 


and using shallow-water theory. The problem of 
oscillation in such a container has been solved by Lamb 
(reference 1, p. 291); the pertinent parts of the solu- 
tion are repeated here in a convenient form. By as- 
suming d = f(p) cos @ (independent of 7), the governing 
differential equation for f can be found from the vari- 


ational principle [Eq. (19)] to be 


p*(p? — 1)f’ + p(3p? — 1)f’ + 
[1 — (2A + 1)p?I|f 


I 


(37) 


or 


[p(p? — 1)f] + [(1/e) — (2A + Iplf = 0 (88) 


Then, appropriate boundary conditions are that 
Note, then, for 
future use, that Eq. (38) implies that the orthogonality 


condition 


1 
f pfm(p)fn(p) = O (39) 
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holds for eigenfunctions f,, and f, associated with dis- 
tinct eigenvalues \,, and i,,. 

The solution of Eq. (38) for the antisymmetrical 
modes is 


n 


fo =D ap*', n= 1,2,... (40) 
j=1 
where 
ay = | 
Qy41 = —([(n? — 7*)/7G + Ia; 
and the corresponding eigenvalues are 
A, = 2n?§ — 1, n=1,2,... (41) 


The first three eigenfunctions, normalized to give 
F121) = 1, ase: 


fi = p 
fo = 3p* — 2p 
fs = 10p®° — 12p* + 3p (42) 


As in the case of the circular canal, the lowest eigen- 
value is A, = 1, corresponding to the frequency of a 
pendulum of length R. It is interesting to compare 
the eigenvalues for the nearly empty canal and spherical 
tank; from Eqs. (41) and (23), 


(2 es _ (eS a + i. 1) 


(3) Sloshing Forces and Motion 


Suppose that the rigid container shown in Fig. 1 is 
subjected to a displacement l(t) in the x direction; 
tank displacements in the y and z direction, as well as 
rotations of the tank about all axes, are assumed to be 
fully constrained. Define a displacement potential Vv 
for the resulting fluid displacement relative to the tank, 
such that the relative displacements are given by 


[OW /Ox, OVW/Oy, OW /Oz} 


and assume W in the form 


¥=> 


a,(U)d,(X, y, 3) 


where the @,’s are the potentials associated with the 
natural modes of liquid oscillation in a stationary, rigid 
tank. Lagrange’s equations will be the basis for 
establishing differential equations relating the general- 
ized coordinates U/(t) and a,(t) to the external force X 
on the tank. The results will be used to calculate the 
“sloshing”’ forces of the liquid on the container wall. 
The potential energy during motion is 


P.E. 


| » f prg(OW/Oy)? dF 
F 


(1/2) f pre(Da.(0¢./dy)]? dF 
F 


where p, is the density. But [see Eq. (5)] 


O¢,/OV = (wn?/gZ)dn 


so that 
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P.E. = (1/2) f pig (d,0,”, ‘g)o,|?> dF (43) 
I 


Note, however, that if ¢,, and @, are modes corre- 
sponding to the distinct frequencies w,, and w,, 


f Vom:VO.dV = 
_ f onV*on,dV + f. onVo-ndS + i) $n (Op,/Oy) dF 
(44) 
Also, by Eqs. (3) and (4), the first two integrals on 


the right-hand side vanish; and so, by Eq. (5), the 
right-hand side becomes 


f (@,,?/2)bmb.d PF 
F 


But then the left-hand side of Eq. (44) is also given by 


(45) 


f (@m?/Z)OmbnidF (46) 
Jf 


and so, for w, # w,, the equality of Eqs. (45) and (46) 


implies that 
J Pm, dF = 0 
- 


f Von: Vo,dV = 0 
. 


(47) 
and also 

(48) 
Hence, the potential energy [Eq. (42)] becomes 


P.E. = (pz/2g) { Xa.t, 1 "dF 
J! 


or, letting a, = { ¢,” dF, 
F 


P.E. = (p2/2g)>5 and, Wn! (49) 
The kinetic energy of the container-liquid system is 


K.E. = (1/2)M.(U)? + (pz/2) X 


J {{U + ¥ a,(0¢,/dx)]? + 
‘ 


[DU an(O¢,/dy)]? + [Qan(Ogn/dz)]?} dV 
where J/, is the container mass. 
Using Eq. (48), together with the energy relation 
for a single mode [see Eqs. (44) and (45) ]. 


(1 2 \vo.|2dV = (o.2/20) f $.2 dF 
V F 


K.E. = (1/2)(M. + Mz) (U)? + 
(pr/2g) Ls anton *(n)* + pn UL) Buds 


f (O¢,,/dx)dV 
mq 


and i, is the total mass of liquid. 


gives 


(50) 
where 


B, = (51) 


ACE 
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The expression for 8,, may be transformed as follows: 


i= fie Ox) [x(0¢,,/Ox)] + 
(0/Oy) [x(0¢,/Oy)] + (0/dz) [x(0d,/0z|] } X 
dV — f. xV*o,d V 
= f x(0¢,,/Ov)dF + f. xVo-ndS 
Finally, Bn = (@,” df xo,dF (52) 


With the Lagrangian L = (K.E. — P.E.) determined 
by Eqs. (50) and (49), Lagrange’s equations 


(d/dt)(OL/OU) — (OL/OU) = X 
(d/dt)(OL/Od,) — (OL /da,) = 0 


yield the differential equations 


a, + w,7a, = —(g/w,?)(B,/a,)U (53) 


(M, + M,)U + ord Bid. = X (54) 


The sloshing force F, of the /iguid acting on the con- 
tainer tn the positive x direction is 
F, = 


S 


MU —X 
Therefore, from Eq. (54), 

F; = —M,U — pr>> Bran (55) 
These results will be cast in convenient forms for the 


circular canal and the spherical tank. 


(3.1) Circular Canal 


The quantities a, and 8, in Eqs. (53)—(55) involve 
only the values of ¢, at the surface, which are propor- 
tional to the eigenfunctions w,, of the integral equation, 


Eq. (14). For convenience, choose w,,(1) = 1; then 


Bn, = (24,2 gar)? f fw(£) dé 
0 


1 
2(aR) f [w(é) 
0 


Introduce, now, the slosh height ¢, at the right-hand 
side of the canal associated with the mth mode. 


|? de 


I| 


f, = a,(d6,/dy) (aR) 
= (d,)(wn?/g) w,(1) 


fag 
= 220n°/ 2 


The total height of slosh at the wall is then >> ¢,. 
Eqs. (52)—(54) may then be modified to 


&, + wn%, = —(\,A)(B,/A,)U (56) 
(M, + Mz)U + 2p,(aR)? >> Bf, = X (57) 
F, = —M,U — 2p,(aR)?>. B,é, (58) 
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TABLE 1 
Eigenvalues for Circular Canal 


r a Aa A;a A3a Ay Ay r 
=i. 6:0 0.0 0.0 0.0 1.0 6.0 15.0 
-—0.8 0.6 0.627 3.23 6.51 1.045 5.38 10.85 
-—(0.6 0.8 0.879 3.98 7.30 1.099 4.97 9.13 
—().4 0.917 1.068 4.34 7.63 1.165 4.7 8.33 
—().2 0.980 1.224 4.56 1.@e 1.249 4.65 7.99 

0.0 1.0 1.360 4.70 7.96 1.360 4.70 7.96 

0.2 0.980 1.482 4.81 8.06 1.513 4.91 8.23 

0.4 0.917 1.596 4.89 8.15 1.742 5.34 8.89 

0.6 0.8 1.706 4.97 8 22 2.13 6.22 10.28 

0.8 0.6 1.822 5.05 8.30 3.04 8.42 13.84 

1.0 0.0 2.018 5.20 8.44 oc o0 G0 
where 

1 
A, = f [w,(&)]* dé (59) 
0 
and 
1 
B, = i tw,(£) dé (60) 
0 


Here, of course, all forces and masses are understood 
to be per unit length of the canal. 

It is important to note that for the circular canal 
these results are valid without the restriction imposed 
in the general development that only rectilinear motion 
of the tank in the x direction was permitted. Rotation 
of the canal about its center produces no sloshing, nor, 
by symmetry, does vertical tank motion produce hori- 
zontal sloshing force. (In fact, it can be shown that 
vertical motion produces no sloshing at all in the circu- 
lar canal.) Since all pressures are in a radial direction, 
the resultant horizontal slosh force F, always passes 
through the center of the canal. ' 


(3.2) Spherical Tank 
In terms of the sloshing natural modes 
o, = v.(r, y) cos kd 


the parameter 8, is 
2aR 2r 
Bn = (wn? »| [ r(r cos @)W,(r, 0) cos k@ dr dé 
0 /7 0 


Thus, it follows that 8, vanishes for all k ¥ 1, justi- 
fying the restriction previously made to consideration 
of only the case k = 1. The values of y, at the free 
surface are proportional to /,(p), where /,(p) = 
(] V p)g.(p),and g,() is an eigenfunction of the integral 
equation, Eq. (28). Normalize /,,(p) so that f,(1) = 1. 
Then, as in the circular canal case, the slosh height 
¢, at the wall is a@,w,?/g, and 


1 2s 
(aR)*(w,? of f p’ f(p) cos® 6dp 
0 0 
1 


m(aR)*(w,,? »| p*f(p) dp 


J0 
1 2r 
a, = (aR)? i plf(p)]? cos? 6 dp 
JoJo 
1 


ma(aR)? f plf(p)|? dp 
0 


B 


II 


Il 


Hence, from Eqs. (53)—(55), 
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&. + w.%, = —(A,a)(D,/C,)U (61) 
(M, + M,)U + rp,(aR)*> Dg, = X (62) 
F, = —M,U — xp,(aR)*> D.E, (63) 
where 
1 1 
C, = f plf.(p)|? dp = [ [g.(p)|?dp (64) 
0 /J0 
and 


1 
D, -f p*t,.(p) dp (65) 
0 


The force F,, again, passes through the center of the 
sphere. 


(4) Numerical Results 


(4.1) Circular Canal 

The first three modes and frequencies for the circular 
canal were found for the values e = 1, 0.8, 0.6, 0.4, 0.2, 
0, —0.2, —0.4, —0.6, and —0O.8, from the 20-by-20 
matrix formulation described in the Appendix. The 
values of (X,,a) and X, (n = 1, 2, 3) thus obtained, to- 
gether with those found analytically fore = —1, are 
given in Table 1. 

The solid curves in Fig. 9 show the variation with e 
of VX (which is proportional to frequency) for the 
first three modes. Worthy of note is the fact that the 
variation with e of the higher frequencies is not mono- 
tonic; the minimum frequency of a given higher mode 
appears to occur slightly below the half full condition, 
whereas the fundamental mode has its lowest fre- 
quency in the nearly empty state. This situation is 
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Fic. 9. Frequencies for circular canal (solid) and spherical tank 
(dashed ). 
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at least partially associated with the fact that, in gen- 
eral, as depth increases, frequencies of free-surface 
oscillation tend to shift:to more closely packed spectra. 
Thus, for example, the ratio \2/A, decreases mono- 
tonically with increasing e. 

For the nearly full case (e = 1), the matrix calculation 
provided the result (Ai@) = 2.018, as compared with 
the upper bound of 2.094 obtained from the one-term 
Galerkin The corresponding asymptotic 
behavior 


solution. 


Vy ~ [(2.018)"/?/(1 — e?)"/4] fore—> 1 


is shown as the dotted curve in Fig. 9; the asymptotic 


behaviors for V2 and WV); are similarly indicated. 
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TABLE 2 
Eigenvalues for Spherical Tank 


€ a Aa Aa A3a Ay Ao PN 
—1 0.0 0.0 0.0 0.0 1.0 7.0 17.0 
0 1.0 1.565 5.34 8.66 1.565 5.34 8.66 
l 0.0 2.78 5.99 9.25 oO x x 


It may be noted that the upper bound A; = 1.367 ob- 


tained by Rayleigh (or by the one-term Galerkin solu- | 


tion of this paper) for the half full canal is very close 
to the value \,; = 1.360 found from the matrix calcula- 
tion. Rayleigh’s improved upper bound, corresponding 
to A; = 1.356, is actually more accurate, but there is 
little doubt that the results obtained from the matrix 
setup are entirely satisfactory for practical application. 
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Mode shapes for circular canal. 
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The modal parameters A, and B, defined by Eqs. 
(58) and (59) are plotted against e in Fig. 10. Except 
for the nearly empty case (e = —1), these parameters 
were found by numerical integration. For e = —1, 
the result [Eq. (22)] for the natural mode shapes 
provides: 

A, 
B,, 


1/(4n — 1) 
i/sforn = | 
= Qforz> 1 


II 


The modal parameter B, is closely related to the 
amount of mth mode induced by lateral acceleration. 
The low values of B, for n> 1 indicate that, in general, 
the higher modes would not have a major influence on 
sloshing forces. 

The mode shapes w,(£) are shown in Fig. 11 for e = 
1,0, and —1. A point of minor theoretical interest in 
connection with the mode shapes is that the slope at 


£ = 1 can be shown to be given by 


w,'(1) = —A,e 


Thus, the slope is negative for e > 0 and positive for 
e < 0; for e = 1, the slope is infinite. The nature of 
these slopes is a direct consequence of the inclination 
of the container wall at the free surface. 


(4.2) Spherical Tank 

The values of (A,@) and X, for the nearly full (e = 1) 
and half full (e = 0) spherical tank were found, for 
n = 1, 2, and 3, from the matrix formulation described 
in the Appendix and are given in Table 2. Also pre- 
sented in this table are the values found analytically 


for the nearly empty tank (e = —1). From the values 
of (A,a) for e = 1, the asymptotic behavior of Vn 
near e = | is shown as the dot-dash curves in Fig. 10, 


and the values of V), for e = 0 and e = 1 are shown 
as the circles in this figure. The dashed curves repre- 
sent estimates of the variation of i, in the intermediate 
ranges, based on the presumably analogous trends for 
the circular canal. 

The modal parameters C, and D, were calculated (by 
numerical integration for e = 1 and e = 0, and analy- 
tically from Eq. (42) for e = —1) and are shown in 
Fig. 12. The curves connecting the calculated points 
in this figure were estimated, with the results of Fig. 9 
for the canal as a guide. 

Finally, the mode shapes f,(p) = (1, V p)en(p) for 
the three calculated cases are given in Fig. 13. 


(5) Comparison With Experiment 


Resonance tests (unpublished) conducted by the 
Lockheed Missile Systems Division on partially filled 
spherical tanks provided the experimental results for 
V«, which are shown as crosses in Fig. 14. The esti- 
mated curve for VA, given in Fig. 10 is reproduced in 
Fig. 14 for comparison. The agreement is generally 
good, and the agreement of experiment with the rigor- 
ously calculated frequency for the half full case is ex- 
cellent. 
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(6) Conclusions 


The results for the circular canal are sufficiently 
complete to be used in the dynamic analysis of such 
canals filled to arbitrary depth. Some degree of cau- 
tion is evidently required, however, in the use of the 
estimated results for the spherical tank in the ranges 
intermediate to the nearly empty, half full, and nearly 
full cases. 


Appendix—Numerical Solution of Integral 
Equations 


The general kernel function [Eq. (13)] for the circular 
canal has a logarithmic singularity at’ = —> 0. For 
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Fic. 12. Modal parameters for spherical tank. 
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f near £, withé + 1, 
G(é, £) = —(1/m) log |— — E| + p(é)  (A-1) 


where 
f 2y(l—€)r Y 
b(t) = — - log- > (A-2 
” 7 "la +7 — (1 — Ors aia 
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Fic. 13. Mode shapes for spherical tank. 
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Near é = 1, 
G(1, —) = —(y/m) log |1 — €| + ylog2. (A-3) ” 


Note, however, that G(0, —) = G(, 0) = 0 and recall 
that w(0) = 0. Then, the integral 


1 
f G(&, &) w(é) dé sc 
0 or 


may be approximated by 


in 
N-1 iA+(A/2) : 7 , . : 
> w(iA) f G(é, &) dE + w(1) [ G(é, &) dé tm 
tak tA —(A/2) J1—-A/2 ne 
(A-4) 
H 
where A = 1/N. Suppose é = 7A; then, ( 
tian =\ JF in 
f G(é, &) dE = AG(jA, 1A) 1# JAN 
t4—(A/2) 
' ae i 
f G(é, &) dé = (A/2)G(j4, 1) JAN 
1—(A/2) 


j4+ (4/2) aa “ 
f G (&, &) dé = Ap(jA) + 


jA—(A/2) HH 
(4/2) 

i] —(1/m) log |t| dt = Aj) p(jA) + 
=(a 


(1/mr)[1 — log (A/2)]} 7 ¥ N 








mé 
% \ oF cal 
| G(1, &) dé = : 
1— (4/2) Vv 
(A/2)} [(y/m) log 2] + (y/m)[1 — log(A/2)|} 
; j ([E 
Hence, the integral equation, Eq. (14), is approxi- 1 
mated by the matrix equation - 
AT }w} = (N/da) }w} (A-5) of 
i spt 
where ;w; is the column vector with elements w( 7A) Th 
(j = 1,2,...N), T is the “integrating” diagonal ma- (% 
trix 
' 2 
1 wh 
1 not 
1/2] 
and A is the symmetrical square matrix with the gen- 
eral element A ;; defined by: 
Ai; = G(tA, jA) ti# IAN 
= p(jA) + (1/m)[1 — log(A/2)} t=y N 
= (y/m)log 2 + (y/m)[1 — log(A/2)] «=jx#N 
Note that, in the absence of a logarithmic singularity, 
the integration scheme would reduce precisely to the 
conventional trapezoidal rule. 
The orthogonality condition for distinct eigenfunc- 
tions w,, and w,, 
{ Wm(E)w,(E)\dE = O 
0 
is approximated, according to trapezoidal rule, by 

















(A-3) 


recall 


&) dé 


(A-4) 


A) + 


gen- 
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‘ity, 
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[wn] T}w.t = 0 (A-6) 


where [w,,] is a row matrix and }w,} is a column matrix. 

The orthogonality condition [Eq. (A-6)] is satisfied 
exactly by the eigenvectors of the matrix equation, 
Eq. (A-5), and, for this reason, the present integration 
scheme is considered to be a reasonably appropriate 
one. 

\n entirely analogous scheme is adaptable to the 
integral equation, Eq. (22), for the spherical tank. In 
the case of the nearly full tank, the behavior of //(p, j) 
near p = pis 


H(p, p’ ~ 
(1/m)log |p — p| + (1/m)(log p + log 8 — 2) (A-7) 
including the case p = I. 
For the half full tank, for p near p, with p + 1, 
H(p, p) ~ —(1/m)logip — p| + 
(1/2) [log p + 2p* K(p?) + log 8 — 2] (A-S) 
On the other hand, /7(1, p) near p = 1 is 
H(1, p) = —(2/m) log (1 — p) + 
(1/r)(2 log S — 2) (A-9) 


Hence, the matrix formulations for the spherical tank 
may be made in a fashion similar to those for the circular 
canal; the 
V pf,(p). 

Eigenvalues and eigenvectors of the matrix C = AT 
{Eq. (A-5)] were found by using a Remington Rand 
1103AF computer on the basis of 20-by-20 matrix 
approximation—i.e., V = 20. The highest eigenvalue 
of C (and, hence, the lowest value of \) and the corre- 
sponding eigenvector were found by matrix iteration. 
from the modified 


eigenvectors are then proportional to 


The second mode was found 


(“swept’’) matrix 
C’ = C — Clu} [w|7 


where w, is the previously determined first eigenvector, 
normalized according to 
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Fic. 14. Experimental frequencies for spherical tank 


After determining the second eigenvalue and eigen- 
vector, again by iteration, C’ was similarly swept to 
provide a matrix C”, direct iteration of which gave the 


third eigenvalue and eigenvector. 


References 


1! Lamb, Sir Horace, Hydrodynamics, 6th Ed., Dover Publica- 
tions, New York, 1945 

2? Graham, E. W., The Forces Produced by Fuel Oscillation in a 
Rectangular Tank, Douglas Aircraft Co., SM-13748, April, 1951 

3’ Lawrence, H. R., Wang, C. T., and Reddy, R. B., Vari- 
ational Solution of Fuel Sloshing Modes, Jet Propulsion, Vol. 28, 
No. 11, November, 1958 

‘ Jahnke, E., and Emde, F., Tables of Functions, 4th Ed., Dover 
Publications, New York, 1945 

5 Grobner, W., and Hofreiter, N., Jntegraltafel, Zweiter Teil, 
Bestimmte Integrale, p. 40, Springer-Verlag, Vienna, 1950 

6 Kellogg, O. D., Foundations of Potential Theory, p. 247, 
Ungar Publishing Co., New York. 

7 Byrd, Paul F., and Friedman, Morris D., Handbook of Elliptic 
Integrals, Springer-Verlag, Berlin, 1954. 





for publication must be a double- or triple-spaced original. 


not acceptable.) Photographs must be on glossy white paper. 
than M/4 in. high. 


and Greek symbols should be clearly distinguished. 





Notice to Contributors 


To expedite review, an additional copy of both the manuscript and figures should be submitted. The manuscript submitted 
There should be wide margins on all sides of the sheets and triple 
spacing around formulas to allow for the marking of directions to the printer. 


Original drawings (jet black India ink on white paper or tracing cloth) should accompany the manuscript. 
The smallest lettering on 8- by 10-in. figures should be no less 


Formulas should be typewritten or carefully written in pen and ink. 


See Inside Back Cover of this issue for further details 


(Blueprints are 


The difference between capital and lower-case letters 

















Viscous Aerodynamic Characteristics in 
Hypersonic Rarefied Gas Flow’ 


RONALD F. PROBSTEIN* ann NELSON H. KEMP** 


Aveo- Everett Research Laboratory 


Summary 


This paper considers the problem of calculating viscous aero- 
dynamic characteristics of blunt bodies at hypersonic speeds and 
at sufficiently high altitudes where the appropriate mean free 
path becomes too large for the use of familiar boundary-layer 
theory but not so large that free molecule concepts apply. 

Results of an order-of-magnitude analysis are presented to 
define the regimes of rarefied gas flow and the limits of con- 
Based on theoretical and experimental evi- 
dence, the complete Navier-Stokes equations are used as a 
model, except ‘‘very close’’ to the free molecule condition. This 
model may not necessarily give the shock wave structure in 
detail but satisfies overall conservation laws and should give a 
reasonably accurate picture of all mean aerodynamic quantities. 

In this ‘“‘intermediate’’ regime there are two fundamental classes 
class and a “‘merged layer” 


tinuum theory. 


of problems: a “‘viscous layer”’ 
class, the latter corresponding to a larger degree of rarefaction. 
For the viscous layer class there is a thin shock wave, but the 
shock layer region between the shock and the body is fully vis- 
cous, although the viscous stresses and conductive heat transfer 
are small at the shock wave boundary. Here, the use of the 
Navier-Stokes equations with outer boundary conditions given 
by the Hugoniot relations is justified. For the merged layer 
class, the shock wave is no longer thin, and the Navier-Stokes 
equations can be used to give a solution which includes the shock 
structure and has free-stream conditions as outer boundary con- 
ditions. A simpler procedure is presented for ‘incipient merged”’ 
conditions where the shock may no longer he considered an in- 
finitesimally thin discontinuity but where it has not thickened 
sufficiently to entail the ‘‘fully merged layer’’ analysis. In this 
case we approximate the shock by a discontinuity obeying con- 
servation laws which include curvature effects, viscous stresses, 
and heat conduction. 

For a sphere and cylinder it is shown that the Navier-Stokes 
equations can be reduced to ordinary differential equations for 
both the viscous and merged layer class of problems. Solutions 
of these equations, when used in connection with hypersonic flow 
problems, are in general only valid in the stagnation region. To 
illustrate the viscous layer solutions, numerical calculations 
have been performed for a sphere and cylinder with the assump- 
tion of constant density in the shock layer, which is a useful 
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approximation at hypersonic speeds. To illustrate the merged 
layer solution, calculations have been carried out for a sphere 
using the incipient merged layer approximation 

Results are presented for detachment distance, surface shear, 
and heat-transfer rate in the stagnation region of a highly cooled 
sphere flying at hypersonic speed. With decreasing Reynolds 
number, the shear and heat transfer are shown to increase above 
the extrapolated boundary-layer values in the viscous layer 
regime and then to begin falling in the incipient merged regime. 
As the Reynolds number decreases in the incipient merged 
regime, the density in the shock layer increases, and the static 
and stagnation enthalpy behind the shock decrease. 

Calculations performed for an insulated sphere show that, with 
decreasing Reynolds number in the incipient merged regime, 
the density in the shock layer decreases; the total enthalpy be- 
hind the shock and at the stagnation point increase so that they 
are higher than the free-stream total enthalpy; and the stag- 
nation-point pressure behaves like the total enthalpy. 

For the highly cooled cylinder in the viscous layer regime, the 
same quantities are presented as for the sphere. The increase 
found in shear and heat transfer above extrapolated boundary- 
layer theory is small, in agreement with vorticity interaction 
theory. 

A discussion is given of the behavior of available experimental 
data for viscous flow quantities in the intermediate regime and 
the behavior predicted by the results of the present calculations 
Qualitative agreement is indicated. 


Symbols 
‘8 = constants in specified surface enthalpy, Eq. (4-13) 

Cp = specific heat at constant pressure 

Cap = component of strain tensor 

ta = boundary-layer constant in shear stress, Eq. (5-20) 

G = vorticity divided by sin 0, Eq. (5-4) 

(gp )Pr=1= boundary-layer constant in heat-transfer rate, Eq 
(5-21) 

h = specific enthalpy 

ho, he = enthalpy functions, Eq. (4-1g) 

hy = bn — hy 

F | = integer specifying number of dimensions, = 0 for 
two-dimensional case, = 1 for axisymmetric 
case 

k = thermal conductivity 

M. = free-stream Mach number 

p = pressure 

po, p2 = pressure functions, Eq. (4-1f) 

rr = Prandtl number, cpu/k 

In» Ig = radial and tangential velocity components 

— qo = stagnation-point surface heat-transfer rate 

R = radial polar coordinate 

R, = radius of curvature of body 


R, = radius of curvature of shock 


Re. = Reynolds number, po UR;/u. 
Re, = Reynolds number, po UR;/jip 

Re, = Reynolds number, po UR, Mo, 
T = absolute temperature 

U = free-stream velocity 

V = Uy-1)Ky4+ 1) 
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= velocity functions, Eqs. (4-1b) and (4-Ic) 

x = distance along body surface measured from stag- 
nation point 

isentropic exponent 

shock layer thickness (detachment or stand-off 
distance of shock from stagnation point) 

6 = shock-wave thickness 

= density ratio across shock, po /p 
polar angle 


d = mean free path 

u = ordinary viscosity coefficient 

u* = dilational or bulk viscosity coefficient 

mM = ordinary viscosity function Eq. (4-1d) 

bo” = dilational viscosity function, Eq. (4-le) 

im = ordinary viscosity function evaluated at an en- 
thalpy of U?/2 

v = kinematic viscosity constant, po /po, 

p = mass density 

p = density function, Eq. (4-la) 

T = shear stress at body 

TaB = components of shear tensor 

’ = dissipation function in energy equation 

y = stream function, Eq. (5-1) 

a = R-dependent part of stream function, Eq. (5-2) 

2 = vorticity interaction parameter, Eq. (5-24) 

Qn = vorticity component normal to flow plane, Eq. 
(5-4) 

Subscripts 

b = conditions at body surface (or wall) 

FM free molecule value 

R,0,n = radial, tangential, and normal directions, respec- 
tively 

s = conditions immediately behind (downstream of) 
shock wave 

stag conditions at stagnation point 

x = conditions in free stream 

Se 3 = differentiation with respect to R 

~ = same order of magnitude (requires dimensional 


consistency ) 


(1) Introduction 


A* INCREASINGLY important problem in aeronautical 
engineering is that of predicting aerodynamic 
characteristics of bodies flying at very high speeds and 
high altitudes. Under these conditions, the flow 
about a given aerodynamic configuration may be suffi- 
ciently rarefied so that the appropriate mean free path 
becomes too large for the use of classical boundary- 
layer theory but not large enough for free molecule con- 
cepts to apply. In such an “intermediate” rarefied gas 
regime, any complete calculation of the flow field 
would in general require the full formulation of kinetic 
theory. At the present time such an approach is not 
satisfactory for obtaining mean aerodynamic character- 
istics such as skin friction and surface heat-transfer 
rate. 

An alternative procedure for making reasonably 
accurate estimates of the mean flow variables at hyper- 
sonic speeds is to postulate the Navier-Stokes equations 
as a model for this intermediate regime. In so doing 
we assume that the flow is not so close to the free 
molecule condition that it is in the first-order molecular 
collision regime, to be defined later. Such a “‘con- 
tinuum’’ approach was suggested by Adams and 
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Probstein,' and the ideas were modified and further 
amplified in Sections 10.1 and 10.2 of Hayes and 
Probstein,? and in reference 3. The assumption that 
the Navier-Stokes model should be a satisfactory one 
for determining mean flow quantities, even when con- 
tinuum concepts may no longer be theoretically justi- 
fied, is empirically based upon certain known experi- 
mental and analytical results which are briefly outlined 
in reference 2. 

Most earlier analytical treatments of the intermediate 
regime have dealt principally with the effect of the 
introduction of slip and temperature jump into the 
This approach to the problem is most adequately 
The present treat- 


flow. 
described in Schaaf and Chambré.‘ 
ment is not meant to imply that slip does not exist as 
an identifiable phenomenon, but only that its effect on 
mean aerodynamic quantities is generally not large 
at hypersonic speeds where our Navier-Stokes model is 
applicable, although, if the body is insulated, it will 
have to be considered in certain cases. 


(2) Rarefied Gas Flow Regimes 


Rarefied gas flow regimes may be classified according 
to the degree of rarefaction as measured by the value 
of the Knudsen number, defined as the ratio of the 
mean free path \ to a characteristic flow length Z. In 
references 1—3, it was emphasized that the concept of 
the Knudsen number must be applied in context, in 
connection with the particular characteristic scale used 
and the particular phenomenon of interest. For ex- 
ample, if we are interested in the flow behind a strong 
shock wave, we must consider the mean free path be- 
hind the shock, whereas if we are studying the structure 
of a strong shock wave, we must consider the character- 
istic length to be its thickness. It is as a result of this 
point of view that we find differences arising in our 
classification of rarefied gas flow regimes when com- 
pared with the original classification of Tsien.° 

We shall summarize, for hypersonic flight speeds, a 
possible classification for the regimes of rarefied gas 
flow. This classification is discussed in greater detail 
in reference 2, Section 10.1. It is based principally 
on consideration of a blunt body, but it is general 
enough so that its main features are applicable to 
slender configurations as well, with appropriately rede- 
fined characteristic lengths. 

In assessing the limits of continuum theory, we shall 
assume that the local mean free path must be small 
compared to the characteristic length of interest for 
continuum theory to be valid. Therefore, we take as 
the basic criterion for continuum flow on an axisym- 
metric blunt body 


rA./A <1 (2-1) 


Here, the subscript s denotes conditions in the shock 
layer behind a shock wave which may not be sharply 
defined, and A is the shock layer thickness (see Fig. 1). 
With R,, the body radius of curvature, the order of A is 
given by (cf. references 1—3) 
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A R, rw € (2 2) 
where 

€ = pu/ Ps (2-3) 


is the density ratio across the shock, taken to be ap- 
proximately constant along the shock and small com- 
pared to one. For air behaving as a perfect gas at 
hypersonic speeds, « — 1/6, while for atmospheric 
hypersonic flight, because of the elevated temperatures 
and resultant activation of more degrees of freedom, 
such as dissociation, ¢ is of the order of 0.1 to 0.07 
or less. At hypersonic speeds the corresponding order 
of magnitude of the pressure p, is p..U* and of the en- 
thalpy, /, is (1/2)U*, where U is the free-stream ve- 
locity. In the shock layer the tangential velocity 
component gy is of the order of Ux/R, with x the dis- 
tance along the surface measured from the stagnation 
point, while the magnitude of the radial velocity com- 
ponent gz varies from el) to U depending on the 
degree of rarefaction. 

Since we are primarily interested in hypersonic flight 
conditions where the temperatures of interest in the 
flow field are high, we may for simplicity assume that 
the viscosity w is proportional to the square root of 
the absolute temperature 7. This corresponds to the 
rigid sphere model of kinetic theory for which Xp is a 
constant. With this assumption, it follows from our 
previous estimates that 


N/A ~d/R, ~ M./Re.« (2-4) 


where the subscript © denotes free-stream conditions, 


M.. is the free-stream Mach number, and Re.. is the 
free-stream Reynolds number 
Re o = Pa UR, M« (2-5) 


We should point out that, with a two-dimensional 
blunt body, these order-of-magnitude estimates would 
have to be somewhat revised since, with ¢« small, the 
shock layer thickness is much larger than on the corre- 
sponding axisymmetric blunt body. Furthermore, 
although as a result of the assumption of hypersonic 
speeds, certain of our free-stream criteria may appear 
similar to those of Tsien,® we point out that these free- 
stream criteria would change with lower speeds and 
different body shapes. 

The preceding order-of-magnitude 
used to define in a rough quantitative manner the 
intermediate regime. We may consider as one limit 
of this intermediate regime the boundary-layer regime 
which exists at Reynolds numbers sufficiently high that 
viscous effects can be taken into account by classical 
boundary-layer theory. The other limit at the very 
low density end is the free molecule regime where inter- 
molecular collisions are neglected and the mean free 
path is everywhere much larger than the characteristic 
length. Based on the classification given in reference 
2, we subdivide the intermediate regime from high to 


very low Reynolds numbers as follows: 


estimates are 


AERO/SPACE 


1960 


MARCH, 


SCIENCES 


(1) Vorticity Interaction Regime 

The vorticity in the inviscid part of the shock layer 
is sufficiently high that boundary-layer concepts must 
be modified. Without external vorticity, the velocity 
at the outer edge of the boundary layer equals the ve- 
locity at the inner edge of the inviscid flow obtained by 
ignoring the boundary layer. With external 
ticity, the velocity at the outer edge of the boundary 
layer must equal the velocity in the inviscid flow at 
the same value of the stream function. A parameter 
characterizing the importance of vorticity interaction 
may be defined as the ratio of the external inviscid 
vorticity to the average vorticity in the boundary 
layer (roughly the velocity at the edge of the boundary 
layer divided by the boundary layer thickness). In 
the case of a blunt axisymmetric body, the inviscid 
velocity profile is approximately linear and the vor- 
ticity interaction parameter is of the order of one. 
From this condition the above order-of-magnitude 
estimates lead to the following rough quantitative 
limit for this regime: 


vor- 


Aa/Ry <e€ (2-6) 


(2) Viscous Layer Regime 


Here the Reynolds numbers are less than required 
for the validity of the vorticity interaction approach; 
the shock layer is a fully viscous continuum amenable 
to treatment with the complete Navier-Stokes equa- 
tions; and the shock wave may be treated as a discon- 
tinuity across which the Hugoniot (shock) relations 
apply. Expressing this condition by requiring the 
shear at the shock (u0q,/O0R), to be small compared 
with the tangential momentum transport at the shock 
(pGe9r)s, We can show that we may use the viscous 
layer concept when 


=<" (2-7) 


(3) Incipient Merged Layer Regime 

Here we consider the shock layer still to be a con- 
tinuum but may no longer consider the shock as a dis- 
eontinuity satisfying the classical Hugoniot relations. 
In one approach to this regime, we determine conditions 
behind the shock by treating it as thin but introducing 
the effects of curvature, viscous stresses, and heat con- 
duction into the integrated conservation equations. 
In this case the density ratio is not given by the usual 
shock-wave relations nor is the tangential velocity 


component constant through the shock wave. This 
regime is delimited by Eq. (2-1), that is, 
Nia/ ly Ke I (2-8) 


(4) Fully Merged Layer Regime 

For this regime the shock layer is taken to be in an 
almost continuum state, and the shock wave is no 
A strict treatment of this regime requires 
We may delimit 


longer thin. 
the full formulation of kinetic theory. 
its range roughly by the condition that 


Aao/Ry < 1 (2-9) 
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(5) Transitional Layer Regime 


The merged shock layer is taken to be in an almost 


first-order collision regime defined below. It is de- 
limited by the inequality 
r . R, > l (2-10) 


(6) First-Order Collision Regime 


In this regime the appropriate local Knudsen number 
is large but not large enough to ensure the full validity 
of the free molecule concept. Here an analysis can be 
carried out by taking into account only collisions be- 
tween a free-stream molecule and a molecule re-emitted 
by the surface. We shall not consider this regime any 
further but instead refer the reader to Lune and Lubon- 
ski,® Baker and Charwat,’?’ Hammerling and Kivel,® 
Liu,’ and Willis.'° 

From our definitions, the flow within the shock layer 
down to and including the incipient merged layer re- 
gime is a continuum and is properly described by the 
compressible Navier-Stokes equations. Although nei- 
ther the fully merged layer regime nor the transitional 
layer regime is a continuum regime, we may for pur- 
poses of obtaining estimates of skin friction, heat 
transfer, and other mean aerodynamic quantities 
postulate the complete Navier-Stokes equations as a 
model for these regimes. This point of view was dis- 
cussed in the Introduction. 

Excluding the first-order collision regime, with the 
Navier-Stokes model there are essentially two classes 
of problems in the whole intermediate regime. In the 
one class we have a discontinuous shock with negligible 
viscous stresses and heat conduction immediately be- 
hind the shock wave, and the Navier-Stokes equations 
are used in the shock layer with outer boundary condi- 
tions given by the Hugoniot shock relations. This 
class of problems we term a “‘viscous layer’’ class, and 
it includes both the viscous layer and vorticity inter- 
action regimes. (See Fig. la.) In the other class of 
problems the Navier-Stokes equations are used with 
outer boundary conditions given by free-stream condi- 
tions. Here, if we choose to include the shock struc- 
ture in the incipient merged layer regime, we can treat 
this regime together with the fully merged layer and 
transitional layer regimes in one class of problems we 
term a “‘merged layer’ class. (See Fig. 1b). We note 
that, if the body is insulated, slip and temperature 
jump will have ‘to be considered in the fully merged 
layer and transitional layer regimes; however, if the 
body is cooled, these factors will not be important in 
affecting the mean flow variables.’ 


(3) Coordinate Systems and Basic Equations 


The approach we shall adopt to solve the viscous 
and merged layer problems allows us, at least for a 
cylinder and a sphere, to reduce the Navier-Stokes 
equations to ordinary differential equations. Solu- 
tions of the resulting equations, when applied to hyper- 
sonic flow past a sphere or cylinder, are in general only 
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applicable to the stagnation region because of the un- 
realistic behavior which is required of such quantities 
as the viscosity, enthalpy, and thermal conductivity 
away from the stagnation region. The procedure for 
obtaining these solutions is to use cylindrical or spherical 
coordinates, respectively, and to seek cylindrically or 
spherically symmetric solutions. Such a procedure is 
analogous to that employed by Lighthill'' in obtaining 
his inviscid constant density spherical solution and that 
used by Whitham!? and Hayes and Probstein® in ob- 
taining the analogous cylindrical solution. These 
treatments made the assumptions of a discontinuous 
shock wave, an inviscid shock layer, a constant density 
ratio across the shock, and constant density in the 
shock layer. Although the approach described in this 
paper is related to the one used by these authors, it is 
more general, particularly for the merged layer class of 
problems. 

An alternative approach to the present one, described 
in references | and 2, is to employ a coordinate system 
of the boundary-layer type and to simplify the Navier- 
Stokes equations, based on an analysis in which the 
shock layer thickness is assumed to be small in com- 
parison to the body radius (or radii) of curvature. 
The actual analysis of the orders of magnitude is identi- 
cal with the classical one for the boundary-layer equa- 
The magnitudes of most of the flow variables 
The results of such 


tions. 
were given in the preceding section. 
an analysis show that the tangential momentum equa- 
tion is the same as the boundary-layer momentum 
equation with a curvature term added. The principal 
difference is that the usual boundary-layer assumption 
that the normal pressure gradient may be neglected 
cannot generally be made, as for example within a 
shock wave. The analysis which is carried out to deter- 
mine the appropriate normal momentum equation is, 
however, the same as in the usual boundary-layer 
theory. Solution of this system of equations in the 
stagnation region of blunt bodies makes use of local 
similarity and takes advantage of the fact that for 
local similarity to apply the shock layer must be thin. 

Both the thin layer approach and the present ap- 
proach employing a complete solution with cylindrical 
or spherical symmetry are in general compatible. As 
just noted, it will be seen that sufficiently far from the 
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stagnation region, that is, for sufficiently large ray 
angles @ (Fig. 2), the complete symmetric solutions are 
in general unjustifiable and may be anomalous. On 
the other hand, for a blunt body the thin layer approach 
is only valid near the stagnation region. A closer ex- 
amination shows that for application to hypersonic 
flight problems the validity of both approaches is 
roughly the same and both require the hypersonic 
condition that the shock layer be thin. However, 
apart from this question, the complete solutions of the 
present paper do afford new solutions of the Navier- 
Stokes equations which may in themselves be of 
interest. 

The basic equations for our problem are the com- 
pressible, heat-conducting Navier-Stokes equations 
written in cylindrical polar or spherical polar coordi- 
nates R and @ with the axis 6 = () directed upstream 
and with respective velocity components gg and qp. 
(See Fig. 2.) These equations for a flow in thermody- 
namic quasi-equilibrium are given in Appendix (A). 

For both the merged and viscous layers with no slip 
or temperature jump, the boundary conditions at the 
body surface R = R, are 


G(R.) = 0, gr(Ro) = 0 (3-la) 


h(Ry) = hy(@) or Oh(R,)/OR = function of @ (3-1b) 


Here h is the specific enthalpy; the subscript b denotes 
conditions at the body surface; and function of @ in 
Eq. (3-1b) is zero for an insulated surface. 

For the merged layer we have free-stream conditions 
as outer boundary conditions: 


(3-2a) 
(3-2b) 


go> Usiné, drm —Ucos0, asR-> @ 


P> Po, P>mpPo, h>h., aRao 


Because of the equation of state, only two of the three 
conditions in Eq. (3-2b) are independent. 

For the viscous layer the outer boundary conditions 
are those existing behind a discontinuous circular or 
spherical shock of radius R,: 


gdr(R;) = —eUcosé@ (3-3a) 


g(R.) = U sin 6, 
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P(R:) = po + poU? cos? (1 — e€) (3-3b) 


A(R;) = ho + (1/2)U? cos? 011 — &) =(3-3e) 
where 


(3-3d) 


o(R,) = p./e 


Eqs. (3-3) are just the Hugoniot relations for an in- 
finitesimally thin shock. We note that ¢ as written 
here is purely a definition, and, to specify the density, 
we must introduce an equation of state between the 
thermodynamic variables p, h, and p. In the present 
paper we shall be concerned with shocks sufficiently 
strong that the density ratio takes on its limiting value, 
which depends only on the thermodynamic state be- 
hind the shock. For purposes of numerical calculation, 
we will use the perfect gas law 


pb = phl{(y — 1)/y] (3-4) 


behind the shock, with y the isentropic exponent there. 
This, combined with the Hugoniot relations, then gives 
for the limiting density ratio 


e= (y7¥ — L/(v + 1) (3-5) 
With Eq. (3-5) it is clear that in this limiting case either 
y or e may be specified interchangeably as a parameter 
of the problem. 

In the viscous layer problem the location of the 
shock wave relative to the body is unknown. Because 
of the free boundary nature of this problem, one addi- 
tional boundary condition will in general be required. 
This point will be discussed further in connection with 
each of the various models of the flow field to be con- 


sidered. 


(4) Symmetric Solutions 


Subject to the boundary conditions stated previously 
we now seek cylindrically and spherically symmetric 
solutions of the complete Navier-Stokes equations 
which are given in Appendix (A). We assume that the 
dependent variables have the functional forms 


p = po(R) (4-1la) 
Jo = u(R) sin 0 (4-1b) 
dr = —v(R) cos 0 (4-1c) 
bh == po(R) cos 6 (4-1d) 
u* = po*(R) cos 0 (4-1e) 
bP = po(R) cos? 6 + po(R) sin? 6 (4-1f) 
h = h)(R) cos? 6 + he(R) sin? 6 (4-1g) 


Here u* is the dilational or bulk coefficient of viscosity. 
For simplicity we replace c,d7 by dh in the energy equa- 
tion and take the Prandtl number Pr = c,u/k to be 
constant. 

With the forms given above substituted into the 
Navier-Stokes equations, we find on equating like 
terms in @ that they reduce to the six simultaneous 
ordinary differential equations 
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(pov) ’ aoe (1 + J) po(u —_ Vv), R = 0 (4-2) 
fo’ + pov’ = 
2[— pov’ + (40* — Ho 3)vpo’ po|’ + 
(1 + j)molu’ — 2v’ — 3(u — v)/RI|/R_ (4-3) 
po’ — pou[(u — v)/R] = 
—(uo/R)\u’ — [(u — v)/R]} (4-4) 
[2(p2 — po)/R] + pot —2u’ + ul[(u — v)/R]} = 
(uo) u’ — [(u — v)/R]f{)’ — 
(4 R)(po’ po)? [uo* — (1 3) po | i 
(uo/R)}(2 + ju’ — (6+ 7)[(u — v)/R]} (4-5) 
— povho’ _ (a, Pr)} (uoho’)’ + (1 + 7) ‘R |uolto’ + 
2[(1 + 7)/R?|uo(he — ho)} — vpo’ + 
(4/3)uo(jv’ + [(u — v)/R]j? — 
(1 — j)v’[(u — v)/R]) — wo*[(vpo’)/po]? (4-6) 
— pithy’ + 2po(u/R)(he — ho) = 
(1 Pr)} (uohts’)’ +- re! + 7) R |poho’ —— 
(4/R2)po(hs — ho)f — vpo’ + 
2(u/R)(p2 — po) + mop u’ — [(u — v)/R]}{? (4-7) 


Here a prime denotes differentiation with respect to 
R, 7 = 0 for two-dimensional flow, and 7 = 1 in axi- 
symmetric flow. We note that Eq. (4-2) derives from 
the continuity equation, Eqs. (4-3)—(4-5) from the 
two momentum equations, and Eqs. (4-6) and (4-7) 
from the energy equation. 

We now have six equations for nine unknowns. 
Two more relations are given by the dependence of the 
viscosity functions on R. For example, these functions 
may be related to the temperature along the stagnation 
streamline 6 = 0. An additional relation which deals 
with the state properties of the fluid is necessary to 
specify the problem. In the illustrative examples 
presented in the following sections, we assume the 
flow in the shock layer to be approximated by one of 
constant density and specify the state properties by 
p = const. The usefulness of this approximation will 
be discussed later. 

If we consider the flow in the shock layer to be one 
of variable density, we must introduce an equation of 
state of the form 


pb = p(p, h) (4-8) 


However, such an equation yields two additional rela- 


tions. For example, if we were to use the perfect gas 
law 

pb = phl{(y — 1)/7] (4-9) 
we would obtain 

Po = poho[(y — 1)/y] (4-10) 
and 

pe = pohe[(y — 1)/7] (4-11) 


These two equations lead to an overdetermined prob- 
lem. The assumed forms of the solution are thus 
seen to be consistent with a kinematic equation of 
State, pp = const., but are not capable of exactly satis- 
fying the added degree of freedom of a thermodynamic 
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relation such as the perfect gas law. We note that, 
in the present case where the stagnation region is of 
principal interest, a useful approximation, for ex- 
ample, may be to drop the thermodynamic function 
h,. This then makes the system determinate, since 
Eqs. (4-7) and (4-11) are also dropped. 

We turn now to the question of whether solutions 
of the form given by Eqs. (4-1) can meet the boundary 
conditions prescribed by (3-1)—(3-3). At the 
body it is clear that the no-slip conditions, Eqs. (3-la), 
If we consider the specification of 


Eqs. 


can be satisfied. 
surface enthalpy, we see Eq. (3-l1b) can satisfy the 
assumed solution forms for either a constant surface 
enthalpy or one which varies as 4 cos? 6 + B sin® 6, 
where 1 and B are constants. We shall not consider 
the specification of 0h/OR at the surface, except for the 
This condition will, 
For the merged 


insulated case where it is zero. 
of course, agree with Eqs. (4-1). 
layer the assumed form of the solutions are compatible 
with the outer boundary conditions, Eqs. (3-2). On 
the other hand, for the viscous layer, our solutions can, 
in general, only be made compatible with the outer 
boundary conditions specified by Eqs. (3-3) when the 
density ratio across the shock can be considered con- 
stant along the shock. For strong shocks at hyper- 
sonic speeds this is a good approximation and will be 
used in the present paper (cf. reference 2). 

We close this section by restating the boundary con- 


ditions in terms of our R-dependent variables. At the 
body Eqs. (3-1) become 

u(R,) = 0, v(R,) = 0 (4-12) 

ho( Rr) = ho, = A, he(R,) =k, = B (4-13) 

ho'(Rr) = 0, he’(Ro) = 0 (4-14) 


where Eq. (4-13) is the boundary condition to be used 
for a specified wall enthalpy, while Eq. (4-14) is to be 
used for an insulated surface. When the constants 
ho, and hy are equal, the surface enthalpy (temperature) 
is constant. 

For the merged layer the outer boundary conditions, 


Kgs. (3-2), become 


u> U, v> U as R— (4-15) 
Pos Po, PomPu, pomp. asR> om (4-16) 
hohe, heh. asR—~> @ (4-17) 


For the viscous layer, the Hugoniot conditions behind 
an infinitesimally thin circular or spherical shock, Eqs. 
(3-3), are 








u(R;) =u, = U, wR.) =vy= «U (4-18) 
Po( Rs) = po, = P« + p0U(1 — €), 
po(Rs) cal pe, = Po (4-19) 


ho(R.) = ho, = ha + (1/2)U%(1 — &), 
h(R) => he, 


h. (4-20) 


Here again we note that ¢ is expressed as a parameter 
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of the problem which is to be specified, and for the 
viscous layer we take it related to y by Eq. (3-5): 


po( Rs) = po, = poe = [(y — 1)/(v+ I) Ie. (4-21) 

As was mentioned earlier, in general, for the viscous 
layer problem one additional boundary condition is 
required. This can be seen by noting that the system 
may be reduced to one of 10th order with ten boundary 
conditions, for example, by eliminating fo) through 
Eq. (4-10). However, because of the unknown loca- 
tion of the shock wave, an 11th boundary condition is 
needed to determine the solution. For this purpose we 
may use the inviscid 6 momentum equation applied 
along the downstream side of the shock, in which the 
Hugoniot values for the velocity and pressure behind 
This yields a specification on 
It corresponds to 


the shock are inserted. 
the first radial derivative of dp. 
specifying the vorticity behind the shock to be that 
calculated from the Hugoniot conditions. We will re- 
fer to the resulting relation as the ‘inviscid’ boundary 
condition. 

In the special case of the constant-density approxi- 
mation, this additional boundary condition is not re- 
quired because the order of the system is reduced by 
one due to the elimination of po’, while the number of 
boundary conditions remains unchanged.* 


(5) Viscous Layer Solution for Sphere and 
Cylinder 


In this section we shall consider the viscous layer on 
a sphere and cylinder. In general, the solution of this 
problem requires the use of the full set of equations 
given in Section (3). In order to illustrate the nature 
of the solutions for the viscous layer, while keeping the 
calculations as simple as possible, we shall restrict our 
considerations to the approximation that the density 
in the viscous layer is constant.t This then approxi- 
mates the case of an insulated body. However, even 
with cooled surfaces, there ts'evicence that the constant- 
density approximation may be a useful one in the stag- 
nation region of blunt bodies at hypersonic flight 
speeds, if we are interested not in the details of the 
flow but in mean surface characteristics such as shear 
stress and heet-transfer rates. An example, for heat- 
transfer rate, is to be found in the comparison of the 
compressible highly cooled stagnation-point boundary- 
layer results of Fay and Riddell'* with the correspond- 
ing constant-density results of, say, Sibulkin.'* This 
comparison shows the constant-density results to agree 
approximately with the compressible ones, provided 


* This fact was first pointed out to the authors by Drs. A. D. 
Wood and G. D. Waldman of the Avco Research and Advanced 
Development Division. 

7 After completion of the present work and the work of refer- 
ence 3, a constant-density viscous layer analysis by Oguchi® 
came to the authors’ attention. His work differs from ours 
principally in that he only considered a sphere and neglected 
His numerical results therefore differ some- 


curvature effects. 
what from those given here. 
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all temperature dependent flow variables in the con- 
stant-density solution are evaluated at an appropriate 
location in the shock layer.”* 

Further justification for considering first the com- 
paratively simple constant-density case can be found 
in the information it provides about the relative changes 
of shear stress and heat-transfer rate as the Reynolds 
number decreases. We may reasonably expect that 
the relative changes in these quantities from their high 
Reynolds number values will not be affected strongly 
by compressibility. t 

As indicated previously, the results obtained are 
principally applicable in a region which js not far re- 
moved from the stagnation point. Therefore, in gen- 
eral, we shall associate the constant density with that 
prevailing on the axis of symmetry behind the normal 
portion of the shock wave, py) since at hypersonic 
speeds in the stagnation region the difference between 
this density and the one at the stagnation point is 
small. With this constant-density 
only the ordinary viscosity coefficient yo enters, since 


approximation, 


the dilational viscosity is everywhere multiplied by 
po’. Again, for simplicity of calculation, we will take 
uo to be constant at its value wo, behind the normal por- 


tion of the shock. Note that u = yo, cos @ still has a 
6 dependence. 

The basic symmetry equations, Eqs. (4-2) to (4-7), 
can be written for this case by dropping the py’ terms. 
The resulting six equations determine the six unknowns 
u, V, Po, p2, ho, ho. We may reduce the unknowns to 
five by introducing a stream function y to satisfy the 
continuity equation. In polar coordinates, y is given 
by 


go(R sin 6)’ = OY/OR, 


gr(R sin 6)? = —(1/R)(Ow/00) (5-1) 


If we write y as 


y = W(R) (sin 6)!+) (5-2) 


and introduce uw and v through Eqs. (4-l1b) and (4-Ic), 


u = W'/R’, v = (1 + j)(W/R'?) (5-3) 


As is usual for a constant-density flow, the momen- 
tum equations can be uncoupled from the energy equa- 
tions by eliminating the pressure and thereby intro- 
ducing the vorticity. In the present case, this is ac- 
complished by differentiating Eq. (4-5) and eliminating 
po’ and p.’ by means of Eqs. (4-3) and (4-4). The re- 
sult is the vorticity equation in the variables u and v. 

The vorticity in the present problem has only one 
component, {,, normal to the plane of the flow and may 


be expressed as 


(1/R)[(ORqo/OR) — (Ogr/20)] = 
}u’ + [(u — v)/R}} sin 0 
(1/R’)[W" + (1 — j)(W’/R) 
(1 + 7)(¥/R?)] sin 6 = G(R) sin 6 


0 
sin 


II 


(5-4) 


II 


t Calculations including the effects of variable density are at 
present being carried out and will be reported on later 
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With this definition of G, the vorticity equation can be 
written in terms of V as 


(R'+4G’)’ — [201 + j)/R!'~](W'/R%)’ + 
(1/v.) [C1 + 7)R°¥(G/R’)’ — (1 — 7)W’G] = 0 (5-5) 
Here, v, = uo,/ po, is the kinematic viscosity evaluated 
on the axis of symmetry behind the normal portion of 
the shock. 
The appropriate energy equations are Eqs. (4-6) and 


f i+7, a wal” ie 3 Pr RW 
i Y he’ | 
ho” + R 1 [ + R S++ OR 


RW 
vy, R? 


A+R 


alan) (4+, 


] 
ho” + a2 hy E re 


ro L(e) ~ (a) (1+ 
2) - 40+ 2a[(gai) J - 


Pr wv 2(1 + 7) - f 
—Ne = —P 
rd a+ 1 r 


The three equations for ¥(R), ho(R), and ho(R) can 
be solved successively rather than simultaneously, 
since one involves only VW and another only VW and hy. 
The vorticity equation, Eq. (5-5), is fourth-order and 
nonlinear in VY, while the energy equations for h. and 
hy are linear and second order. In fact, Eq. (5-7) can 
be treated as a linear first-order equation in fy’. Once 
WV and its derivatives are found by solving Eq. (5-5), 
po(R) and po(R) — p2(R) can be obtained from the con- 
stant-density (4-3) and (4-5), the 
former requiring a single quadrature on R. 

From Eq. (4-12) the boundary conditions for the 
vorticity equation at the body are 


versions of Eqs. 


V(R,) = 0, W'(R,) = 0 (5-8) 
while at the shock, from Eq. (4-18), 
V(R,) =vV, = 
[leU/(1 +7) |R.'%, W'(R,) =v.’ = UR, (5-9) 


One more boundary condition in addition to these 
four is needed for the solution of the differential equa- 
tion, Eq. (5-5), because of the unknown location of the 
boundary at R. However, it should be noted that this 
does not represent an additional boundary condition for 
the basic constant-density differential system, such as 
is required in the variable-density model. This rela- 
tion must be introduced here only because the pressure 
functions pp and pf. are eliminated by cross differenti- 
ation, which reduces the order of the differential system 
by one while resulting in the loss of the two pressure 
boundary conditions at the shock. This necessitates 
the introduction of one additional relation. We ob- 
tain the appropriate relation by satisfying the 6 momem- 
tum equation, Eq. (4-5), along the shock, where V(R,) 
and W’(R,) are given by Eq. (5-9) and where 2(p2. — 


2 V 


(4-7) with po’ = 0. A discussion of the form of the 
energy equation to be used with the assumption of con- 
stant density can be found in reference 21. The pres- 
sure terms may be eliminated by using the constant- 
density forms of Eqs. (4-3)—(4-5). Inspection of the 
constant-density energy equations shows that 
However, 


two 
they are coupled. 
the other and the new enthalpy variable he = he — Io 
is introduced, the resulting equation depends only on 
ho and WV. It, and the constant-density version of 
4-6), become, using the stream function WV: 


)}+ 


if one is subtracted from 


Eq. ( 


]- -m le) -(e 


. Wy’ 1+) lw 
) + Ri (ei) + R G Ot >, za si 
 - vw \’ +jvV w\) 
idea igen) + | zw) +(e) I ar 
4(1 + 27 | ( = )} + 
J Rit 
., v v \’ l+jwv w\) .- 
l ~4 
( + J) R?2+) | (ax) (: . Vv. =) ¥ (F) i} o 


po) = (0p/08)/(sin 6 cos 6) is obtained from the Hugo- 
niot relation Eq. (3-3b). The resulting relation be- 
tween V’’’(R,) and ¥"(R,) is 


V,’"R2-/(U Re,) = 
—{{(1 — j)/Re.] + 1}(¥."R,!-/U) + 


[7 — 3 + 2e + (1/e)] + [4(1 — ©)/Re,] (5-10) 


We have used the subscript s here to indicate conditions 
behind the shock. The Reynolds number Re, is de- 
fined by 


Re, = p. UR, Mo, (5-11) 


and is appropriate for the viscous layer problem as 
treated here. We shall term Eq. (5-10) the ‘‘viscous” 
boundary condition for Eq. (5-5). 

For a constant surface enthalpy the wall boundary 
conditions on the enthalpy equations, from Eq. (4-13), 
are 
(5-12) 


ho(R,) = ho,, ho(R>) (he — ho), = 0 


With the hypersonic strong shock approximation (MM ..” 
cos? 6 > 1) we may neglect /.. in comparison with h, 
and we may neglect e®. From Eq. (4-20) this leads to 


the boundary conditions at the shock 
2)U?, he(R.) = (ho — he). = 


ho( Rs) = (1, 
(5-13) 


Of course, even at hypersonic speeds for any finite 
Mach number with sufficiently large @, the strong shock 
approximations will break down. However, as al- 
ready pointed out, our considerations are in general 
restricted to the stagnation region where the approxi- 
mation is valid. 
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For numerical integration, Eqs. (5-5)—(5-7) may be 
made nondimensional by using U as the unit of velocity, 
R, as the unit of length, and (i/2)U? as the unit of 
enthalpy. Then, the dimensionless parameters of the 
vorticity equation and boundary conditions are e« and 
Re,, while the enthalpy equations and appropriate 
boundary conditions introduce Pr and /ho,/(1/2)U? 
for the constant wall enthalpy case. Once e and Re, 
are chosen, the vorticity equation is integrated as a 
two-point boundary value problem starting at the 
shock. The stream function variable WV, and its 
derivative VY,’ are known from Eq. (5-9). A value of 
V,” is assumed, and V,’’’ is computed from Eq. (5-10). 
The integration then proceeds toward the body, whose 
initially unknown location is to be determined by the 
simultaneous satisfaction of the no-slip conditions at 
the sarne value of R. This will occur for only one 


value of ¥,”, so that the integration must be repeated 
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with different assumed values of VY,” until the correct 
value is found, thus determining R,. Once the vor- 
ticity equation is solved, the enthalpy equation for 
hy/(1/2)U? is integrated starting from the value —1 
at the shock. Various values of ho’ are assumed until 
ho, = ( is satisfied at the now known body location. 
Finally, the enthalpy equation for /o/(1/2) lL” is solved 
in the same way, using the value | on the shock and the 
prescribed value of /,/(1/2) U* on the body. 

The procedure described has been carried out on an 
IBM 650 digital computer for both the sphere and 
cylinder cases. The parameters were chosen to corre- 
spond to a highly cooled body (surface temperature of 
the order of the free-stream temperature) in atmospheric 


hypersonic flight (kh, ~ U?/2). In particular, we took 


e= 0.1, Pr = 0.71, M,/(1/2)U? > 0 (5-14a) 


For large wall cooling the exact value of /y,/(1/2)U* 
is not important as long as it is sufficiently small. As 
we pointed out before in Eq. (3-5), for the case of a 
strong shock wave in a perfect gas where the density 
ratio approaches its limiting value, we may write this 
limiting density ratio in terms of the isentropic ex- 
ponent behind the shock—1.e., « = (vy — 1)/(y + 1). 
For e = 0.1, the corresponding value of y is then 


(5-14b) 


The use of y as the parameter rather than ¢ is somewhat 
more appropriate, since, for the incipient merged layer 
to be treated in the following section, ¢ is not inde- 
pendent of Reynolds number and we must specify the 
gas through y. 

The range of Reynolds numbers Re, for the integra- 
tion was taken to be between 10 and 1,000 with the 
upper limit chosen at a point where the viscous layer 
results could be expected to agree closely with the re- 
sults of boundary-layer theory. The lower limit was 
chosen somewhat below the Reynolds number for 
which the viscous layer theory could be considered 
valid. According to our estimate of Eq. (2-7), the 
viscous layer analysis is valid for axisymmetric flow 
so long as 


(eRe,)1 < 1 (5-15) 


In arriving at this form of the inequality, we have used 
the order-of-magnitude estimates 

Meo /Me ™ (T /T;)' 2~ j?MU. 
For a density ratio e = 0.1, Eq. (5-15) implies that the 
viscous layer model breaks down at a Reynolds number 
somewhat less than 100. 

Fig. 3 shows the calculated shock layer thickness for 
the sphere and cylinder. For the sphere two values 
are indicated: one for the “viscous” boundary condi- 
tion specified by Eq. (5-10) which is correct for our 
model, and the other for an “‘inviscid’’ boundary con- 
dition which is the specification of V’(R,) through Eq. 
(5-10) by allowing Re, > ©. This “inviscid” condi- 
tion corresponds to that suggested at the end of Section 
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(4) to be used as the additional boundary condition in 
the variable-density case. The use of such an “‘in- 
viscid’ condition in the constant-density case, although 
not self-consistent, can in part serve to indicate when 
the viscous layer model itself is breaking down. This 
follows from the fact that our model is only applicable 
when the viscous shear and heat conduction immedi- 
ately behind the shock are negligible. Any appreciable 
differences in the results found by using the two differ- 
ent conditions would indicate that the assumptions 
of our model are not being met. It can be seen that 
divergence between the results begins to occur at a 
Reynolds number somewhat below 100, which is the 
approximate value our order-of-magnitude analysis 
would indicate as a lower bound for the validity of the 
axisymmetric viscous layer analysis. Similar results 
are found for the other flow variables. These calcula- 
tions thus bear out our remarks concerning the use of 
the viscous or inviscid shock boundary condition. 
All results presented in the remainder of the paper are 
for the correct viscous boundary condition. We again 
caution the reader that the order-of-magnitude esti- 
mates must be altered somewhat for the two-dimen- 
sional case since A/R, is no longer of order «. Never- 
theless, the range of validity of the analysis for the 
cylinder will be about the same as for the sphere for 
most of the flow variables. 

From Fig. 3 for both the sphere and cylinder, the 
shock layer thickness may be seen to increase with de- 
creasing Reynolds number. With the present constant- 
density approximation, this increase results from the 


increase in the usual displacement thickness, which is 
included in A in a viscous layer analysis, plus the added 
displacement effect resulting from the inclusion of vis- 
cous terms in the shock layer as a whole. (If com- 
pressibility with large wall cooling had been taken in- 
to account, the ordinary displacement thickness would 
be negative and so tend to reduce A below its inviscid 
value at high Reynolds numbers where usual boundarv- 
layer theory is valid. However, the inclusion of the 
viscous terms in the whole shock layer, and not just 
near the wall, would still tend to make A increase with 
decreasing Reynolds numbers even in the compressible 
case.) On Fig. 3 we have also shown for reference the 
infinite Reynolds number stand-off distance correspond- 
ing to the calculated result for the Lighthill'! inviscid 
sphere solution and the Whitham’? and Hayes-Prob- 
stein’ inviscid cylinder solution. 

On Figs. 4 and 5 are plotted, respectively, the surface 
shear and heat-transfer rate results for the sphere and 
cylinder. Both the shear and heat-transfer rate are 
made dimensionless with respect to the free molecule 
value for infinite Mach number and completely diffuse 
reflection (all accommodation coefficients unity). These 
values are 

Trm = p.U* sin 6 cos 6 (5-16) 


— Gru = (1/2)p.U* cos 0 (5-17) 


The viscous layer curves for the shear represent the 
shear at any point on the body, since 7 has only a 
component proportional to sin @ cos 6, that is 
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TABLE 1 
Inviscid Flow and Boundary-Layer Constants 
——Constant-Density, « = 0.1 = 
——R./Ri—— _ U-"(dqg/d8 stag 
Approxi- Com- Approxi- Com- 


mate plete mate plete i an (0) ee 
Cylinder 
(j=0) 1.15 1.17 0.548 0.576 1.233 0.570 
Sphere 
(j=1) 1.11 1.07 0.516 0.518 0.928 0.540 
t = [u(0ge/OR) |, = wo,u,’ sin @cos@ (5-18) 


On the other hand, the heat transfer is given by the rela- 
tion 


—qg = [Rk(OT/OR)), = 
(uo,/Pr) [ho,’ + (he — ho)» sin® 0] cos 6 = 


— [do + ge sin? @] cos @ (5-19) 


and we have plotted only the stagnation-point com- 
ponent — qp. 

For comparison purposes we have also shown the free 
molecule values from Eqs. (5-16) and (5-17) and 
the constant-density boundary-layer stagnation-point 
values. These stagnation-point boundary-layer results 
may be expressed as 


tRe./trm = [(1/U)(dqe/d8) Jstag(Re R,)'” x 
(Re,/«€)'/7(1 + jf” 


—GoRe./—qrau = [(1/U)(dqo/dé) Mag(Re RR)" x 
(Re./e)'/?[(1 + j)'?/Pr]Pr4(gp") pr =1 


(5-20) 


(5-21) 


where /,,” and (g,’) p,-; are constants obtained from the 
appropriate boundary-layer solutions. The inviscid 
values of the stagnation-point velocity gradient (dg, + 
0) ing and R,/R, can be found from the complete in- 
viscid constant-density solutions’."!,!" or can be cal- 
culated approximately from the relations (see Chap. 4, 
reference 2) 


[(1/U) (dqo/d9) leae = {[(9 — j)/B]e}1 (5-22) 


R,/R. = 1 — e{(1/2) In (4/36) ]!~? (5-23) 


The approximate relations are given in order to show 
roughly the « dependence. In Table 1 are presented 
the values of the velocity gradient and (R,/R,) for « = 
0.1, as well as the appropriate boundary-layer constants 
(see, for example, references 14 and 15). 

An interesting feature of the viscous layer results 
shown in Figs. 4 and 5 is that for a given Reynolds 
number the shear and heat transfer are higher than 
would be given by an extrapolation of boundary-layer 
theory to low Reynolds numbers. This result is con- 
trary to previous conjectures by many people, including 
one of the present authors (see Fig. 2 of reference 1). 
Shear and heat-transfer values higher than extrapolated 
boundary-layer theory cannot persist at indefinitely 
low Reynolds numbers, since they will eventually be- 
come higher than free molecule values. Adams and 


Probstein' have shown that using the Navier-Stokes 
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equations it is not physically possible on a cold blunt 
body to exceed the free molecule shear stress and heat- 
transfer rate. The fact that the viscous layer results 
on Figs. 4 and 5 do become higher than the free mole- 
cule values near Re, = 10 can be attributed to a break- 
down in the viscous layer model at sufficiently low 
Reynolds numbers. This model does not include 
curvature, viscous, or heat conductive terms in the 
conservation equations used across the shock wave. 
As the Reynolds number decreases these terms within 
the viscous layer make increasingly larger contribu- 
tions to the flux of momentum and energy toward the 
body. We find asa result, because of the discontinuous 
shock-wave model, that more momentum and energy 
flow out of the shock toward the body than into the 
shock from the free stream. This leads eventually to 
an increase of shear and heat transfer above the free 
molecule values. Ata Reynolds number Re, somewhat 
below 100, where curvature, viscous, and conductive 
terms in the shock conservation laws become impor- 
tant, we expect the viscous layer model to become in- 
valid and a merged layer analysis to become necessary. 
We may anticipate that such an analysis would provide 
a drop in the shear and heat-transfer curves relative 
to the viscous layer values. The results of an analysis 
which do indicate this drop are presented in the follow- 
ing section. 

The fact that the first departure from boundary- 
layer theory is toward higher shear and heat transfer, 
as shown by the viscous layer results in Figs. 4 and 5, 
has been indicated previously by vorticity interaction 
theory (see reference 2, Section (9.6)), where the outer 
boundary condition of the momentum equation is 
modified to account for the strong vorticity present in 
the inviscid field. Such an analysis is appropriate to 
the vorticity interaction regime defined in Section (2). 

On Figs. 4(a) and 5(a) we have shown for comparison 
some incompressible, axisymmetric, stagnation-point 
vorticity interaction results. These results were ob- 
tained by using some calculations by Kemp”! in which 
the classical constant-density boundary-layer momen- 
tum and energy equations were integrated with the 
vorticity interaction boundary condition (/’)? ~ 1 + 
20f on the momentum equation instead of the usual 
condition f’ — 1 (ef. reference 2, p. 371). Here f is the 
usual boundary layer reduced stream function and & 
is the vorticity interaction parameter. The results 
of these integrations for 2 < 1 may be expressed ap- 
proximately as follows: the usual shear function f,” is 
multiplied by 1 + 0.782 and the usual heat-transfer 
function g,’ is multiplied by 1 + 0.360. 

For the constant-density case we are considering, & 
is proportional to the constant vorticity in the inviscid 
flow divided by the average vorticity in the boundary 
layer and may be written as 


@ = [(1/U)(dgo/d8) Jai2z” ((Ro/U) (0"qo/OR2) |stag X 
(¢/Re,)'/"(R,/2Ry)"”? (5-24) 


The quantities with subscript stag, along with R,/R,, 
are evaluated using the inviscid constant-density sphere 
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solution of Lighthill.'! All of these quantities have 
been given previously for « = (0.1 in Table 1, except 


[(Rp/U)(07qo/OROP) |stag = 
(R,/R:)?((1 — ©)?/e] — [(1/U) (dgqo/d8) | stag 


(5-25) 


which is equal to 6.035 for « = 0.1. Using these values 
we find 2 = 3.75/(Re,)'”*, which was the result used to 
obtain the points on Figs. 4(a) and 5(a). 

The good agreement between vorticity interaction 
theory and viscous layer theory down to Re, ~ 300 in 
the axisymmetric case indicates that it is primarily the 
strong vorticity in the essentially inviscid part of the 
flow behind the shock wave which accounts for the 
initial departures from boundary-layer theory. The 
viscous layer theory, which includes terms in the 
Navier-Stokes equations not included in vorticity inter- 
action theory, carries us down about another factor of 
3(~ e '/*) in Reynolds number, in agreement with the 
order-of-magnitude estimates of Eqs. (2-6) and (2-7). 
In this regime, the shear and heat transfer are some- 
what higher than vorticity interaction predicts. 

The agreement of viscous layer and vorticity inter- 
action results for the sphere at the higher Reynolds 
numbers also leads to a better understanding of the 
cylinder calculations where the viscous layer effect is 
seen to be smaller than for axisymmetric flow at a given 
Re,. This is perfectly consistent, since it can be shown 
that for the cylinder, vorticity interaction is a second- 
order effect (in Reynolds number) as compared to a 
first-order effect for the sphere. We may put this 
another way by noting that at the surface of the cylin- 
der, for the inviscid constant-density case, the vor- 
ticity is zero. For a sphere under the same conditions 
the surface vorticity has a relatively large finite value 
since the velocity profile in the inviscid layer is approxi- 
mately linear in the radial coordinate. Therefore, 
equal departures from boundary-layer theory occur at 
a much lower Reynolds number for a cylinder than for a 
sphere. 

A further confirmation of the ideas we have ex- 
pressed is found in the calculations of the pressure dis- 
tribution through the shock layer, which show that for 
Reynolds numbers greater than 300 the layer is essen- 
tially a constant-pressure layer in the neighborhood 
of the stagnation region. This is immediately recog- 
nized as one of the basic boundary-layer assumptions. 
At these higher Reynolds numbers the calculated be- 
havior of the velocity profiles for both the cylinder 
and the sphere is also found to be consistent with the 
vorticity interaction picture just given. 


(6) Incipient Merged Layer 


In the preceding section we saw that for a Reynolds 
number Re, somewhat below 100 the viscous layer model 
is no longer satisfactory and we may no longer consider 
the shock as a discontinuity following the classical 
Hugoniot relations. We can now proceed to solve the 
fully merged layer problem, including the shock struc- 


ture. Alternatively, we may consider the shock slightly 


thickened, though still essentially a discontinuity, but 
with curvature effects, viscous stresses, and heat con- 
duction taken into account in a manner which satisfies 
overall energy and momentum conservation laws 
through the shock wave. Because it is the simpler 
of the two approaches, we choose the latter course, al- 
though we recognize according to our estimate of Eq. 
(2-8) that such an incipient merged layer solution can 
at best carry us a factor of three to five lower in Reyn- 
olds number. The importance of such a solution, 
however, is in its alteration of the behavior of the vis- 
cous layer solutions so as to bring the shear and heat 
transfer down in the manner indicated in the preceding 
section. 

Basically this procedure consists of modifying the 
shock boundary conditions for the viscous layer solu- 
tion. Determination of these modified boundary 
conditions involves obtaining integrated conservation 
laws for a shock wave which include viscous and heat 
conduction terms as well as curvature terms.* Our 
method consists of integrating the Navier-Stokes 
equations from behind the shock, which to some extent 
is arbitrarily defined, to upstream infinity where it is 
assumed that all flow gradients vanish. In order to 
simplify the equations, the orders of magnitude of the 
remaining terms are estimated using values consistent 
with the Hugoniot equations, and terms of higher order 
are dropped. The remaining integrals are then evalu- 
ated using the one-dimensional shock-wave solution 
(including viscosity and heat conduction). With these 
modified boundary conditions the density ratio can no 
longer be specified in advance, nor can the tangential 
velocity component be considered a constant through 
the shock. 

The basic form of the Navier-Stokes equations to be 
used has already been given in Eqs. (4-2)—(4-6). These 
equations are repeated in Appendix (B), Eqs. (B-1)- 
(B-5), in a form more suitable for integration. The 
dilational viscosity coefficient is taken zero for sim- 
plicity, although it could be included in the analysis. 
These equations are integrated with respect to R from 
R, to ~, and free-stream boundary conditions along 
with the additional criterion that all gradients vanish 
are then applied at upstream infinity. The resulting 
integral relations are given in Appendix (B), Eqs. (B-6)- 
(B-10). These equations, together with the gas law, 
replace the Hugoniot relations and are used as the 
boundary conditions behind the shock for the following 
quantities: v, €, po, po, u, u’, ho, ho’, he. For the 
moment we note that Eq. (4-7) for 42 and its derivative 
is not considered as it would lead to the boundary con- 


* After completion of the present work it was brought to the 
authors’ attention that the idea of modifying the Hugoniot rela- 
tions in the stagnation region behind a detached shock wave to 
take into account curvature, viscosity, and heat conduction 
effects in the flow behind the shock was proposed in reference 26. 
However, in that work terms dependent on the thickness of the 
shock wave were neglected, and, since these terms are of the same 
order as those taken into account, the theory presented is in- 
consistent; this casts doubt on the results. 
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ditions for the flow behind the shock wave being over- tegrands all vanish as R— ©. The range of R over 
specified. This is analogous to the overdetermined which the integrands have a significant value is taken 
system discussed in Section (4) in connection with the as the shock thickness 6,. In the integrands, quanti- 
compressible solutions. ties are estimated by using the largest value they take 

The integrated Eqs. (B-6)—(B-10) contain terms to on. For all quantities but v, which is taken of order 
be evaluated behind the shock wave at R,, and integrals U, this leads to the same estimates as were employed 
through the shock wave. The orders of the former above for the integrated terms. Since v and fy change 
terms are estimated using the Hugoniot conditions, rapidly through the shock wave, their derivatives can 
Eqs. (4-18)—(4-20), with p. neglected in comparison be large, and for them the operator d/dR must be taken 
with p.U* and h.. neglected in comparison to U?/2. as 6,~'. On the other hand, because u does not change 
In addition we take uo ~ wo,, R ~ R,. Where deriv- rapidly, wu’ may be taken the same as for the integrated 


terms. After all terms in Eqs. (B-6)—(B-10) have been 
estimated in this manner, we drop terms which are 


smaller than Re,~! or 6,/R, times the largest terms in 
each equation, remembering that ¢ is a small quantity. 


atives appear in these terms the operator d/dR is esti- 
mated to be of the order of magnitude of A~!, with A 
of the order of «R, according to Eq. (2-2). In the in- 
tegrals in Eqs. (B-6) to (B-10), we note that the in- 


With « = p./py,, the resulting conservation equations 
are: 
Us - ” Pov Uu dR 
mera f(t 1) 2] - 
l ‘| . I. pol \v R 
; UP ; - pv u—vdR 41 /R,u,’ l | Us 
. io — -— 0 +a = ese = ( ~ + = ) (6-2) 
po? «U? Rs Poll l R 3 Ke, l 2 
25 - pou (u 1R 
- Sue f po  ( ~~ )' (6-3) 
pat R, PoU? \v R 


8s 


Us ( hos 4 “) - 1 1 ; Rho,’ 1 op j fo Ho Ri’ dR 
eU 2p Re, (1/2)U*Pr ~ Re,Pr Jr, wo, (1/2)U?2 R 
Rg 


(1/2) U? UU? 
4(] + j) {- Mo uv ov dR 4 4 l + J f. Ho R,v' u—vdR 
Re. Jr, mw,U U R ' 3| Ra Jr, mw, U U R 


Gl +7  — ( ho 4 “\(° .) dR *% 
( wu -) 
Af, poU \(1/2)U? U?/\v R : 


In these equations the terms remaining are exact as written. We now proceed to approximate them in a manner 
consistent with the order we are considering. In the integrated terms which are multiplied by Re,~', we need only 
evaluate quantities to zero order, since terms of order higher than Re,~! have already been dropped. Therefore, we 
put uw, = U from the inviscid Hugoniot relation, and Ry,’/U = (1 + Jj) from the zero order inviscid constant-den- 
sity continuity equation [Eq. (4.2) with p. = constant and the v term dropped]. Of course, since we are using these 
equations to obtain boundary conditions on u,’ and /o,’, these quantities are left undetermined. We note that v,/eU 
which appears in Eqs. (6-4) and (6-5) is to be evaluated using the continuity equation (6-1), since it multiplies a term 
of orderone. On the other hand, in Eq. (6-2), v./eU/ multiplies v, / U, which is of order ¢, so there we must set v,/¢U = 1. 

We now approximate the integrals to consistent order by employing the hypersonic one-dimensional shock-wave 
structure equations, which can, for example, be obtained from the above equations in the limit R—~> ©. In par- 


Uglls ~ £2 pou (u dR ° po dR inn” 4 ] = uo Rv’ dR : 
pote Fe th. ee -aclt e - 3) -: f << (6-4) 
eV? Re Pol? \p R on Ss oe he: l 3 Re, Jr, mo, UR 


ticular, we use the results 
Pos¥v, = rm oe Po + pov” = Pa U? — (4 3) wor’ 


along with the fact that wu = U. In addition we take R = R,, since R only changes through the shock by 6,._ The 


resulting equations are 


v,/U = (1 — (1 +A [(U/v) — 1](dR/R,) t (6-6) 
Rs 
bo, 2, re ( 4 dR Al+ s)  2AL+ 9) 7’ mo do PAs 
#1-—-(i Ss ee eine = 6- 
pa? ," ere I. U/ R, a. ” -/- I me U = 


: . ny ag (© - 1) = a (1 *) 6-8) 
paU? J. v R, 1+ fj eU a 
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Us wt i dR ° v\ dR 8 1 U yo dv 
aif -1)F+2f oes eet | —_—- 
Rs \U R, Rs U/ R, 3 Re, Jy pmo, U 
3 (Ea. 
Re, U 


R, ho, , 


ho, oe" 1 


1/2)02 * U? ~ Re, (1/ 2) U2Pr 


II 


“2 5 Ho _ ho _ 4 

+ Re,Pr 2 mo, (1/2)U? © 
145 (° m (1-2) 24 td 
3 Re, V Mos U/ U Re, 


(+; ( Nos 4 “) ss (S 1) dR 
Ni\agjwet v I. " R, 


Here the lower limit V = U(y — 1)/(y + 1), the Hugoniot value of 2, 
To evaluate the remaining integrals we must first specify certain fluid properties. 
which leads to Becker’s energy integral for a one-dimensional shock, 


ho + (1/2)v? = const. 


In addition we assume the viscosity uo to be proportional to fy! 
model. 
velocity profile. (See, for example, Hayes," p. 456.) 
imation and drop terms leading to higher order contributions in Re,. 
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(6-9) 


t wo dv 
v wo, U 


U Mo UV dv 
-— of 
f, wo, U L 
1) =| (6-10 
a 


In particular, we take Pr = 3/4, 


(6-11) 


2, which is consistent with our perfect gas rigid-sphere 
With these fluid property assumptions, we are led to a first-order differential equation for the shock-wave 
In writing this result we again make the hypersonic approx- 
The resulting differential equation is 


(6-12) 


It is well known that this velocity profile leads to a shock thickness which is logarithmically infinite and so is not 


suitable for use in approximating the integrals in which v appears in the integrands. 


We therefore follow the usual 


practice for defining a shock thickness and use a linear velocity profile whose slope is the maximum value of dv/dR 


With (y — 1)/(y7 + 1) small, this slope is 


(=) _3 vy U/R, ‘ft y (7 iy] 
mM as41 Pe “\y +1 


The velocity profile to be used is then written 


oe (dv/dR) maz(R axed R,) + z 


obtained from Eq. (6-12). 


From this profile, the shock thickness is defined by setting vy = U at R — R, = 4,, which yields 
6. = — 25)/ (do, 


Using these equations, the following integrations can be performed: 


. (f )s U/R, ( U ) 
f —-—]1 —_- = In- — 1 
Ry \U R, (dv/dR bi Vs 
oe (: “) dR_ 1 _ U/R, 
J. U] R, 2 (dv/dR) maz 


Here we have dropped all terms of higher order than Re,— 
The remaining integrals are evaluated using Becker’s integral Eq. (6-11) and the relation jo 


are 
Up vdv 1 f¥ (: _ dv 
I. Ho, UU 4 i U? U? 3 
a ?.. 
v Mo, U 
[mt ff (om dle 
caus ma, (1/2) U? aya: \(1/2)U?7 (1/2) U® 3 


where again we have dropped terms of higher order. 
Using the preceding results, the modified boundary conditions can be written 


GR) maz 


(6-13) 


(6-14) 


(6-15) 


(6-16) 


(6-17) 


~ ho. The results 


(6-18) 


(6-19) 


(6-20) 
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Us J 8 ¥ l mal (7 =e am U ) 
= eft — 5 Af +2(2—) fin - - 1) - 
- ( 3y+1 Re, + y+1 = m f (6-21) 
Po, Vs 2(1 + 7) J ia 8 7 | (7 = "yt 
' ieee ia a 1+2 fais 
p.U? l Re, l 4 af 3y+1 + y+1 { (6-22) 
po, (1 *) 
p.U2 1+ ” (6-23) 
p_U? '+3 el 
7] l Vs 1 J Ss Y | (2 ame ?) 1 R,u,’) 
eee = l- 7+ 1+2 “ ba 
l ey x) +g \ a syt+1 + +) us (6-24) 
No, Us" 1 F oa 4 Rho, , ) 
vig at 3 (1 + J) i. ; (6-25) 


=—_ l —_ « 
(1/2) U? U' 


These equations are appropriate shock boundary conditions for the incipient merged layer. 


3. 3q/2)Us 


They replace the 


Hugoniot relations with p.. and /.. negligible and reduce to them for Re, ©. They are not all independent since 
Po( Rs), po( Rs), and io(R,) are assumed related through the perfect gas state equation, which from Eq. (4-10) may be 


written 


ho,/ (1/2) VU? = [2y/(v — 1) Je(Po0,/p @U?) (6-26) 


Similarly from Eq. (4-11), we can express the corresponding condition on h.(R,)._ We write this condition in terms 


of our dependent variable /2 as 


he,/(1/2)U? = [2y/(y — 1) ]e[(b2, — ho,)/p 2U?] (6-27) 


To integrate the vorticity equation, an additional condition at R,, corresponding to Eq. (5-10), is required. 


This is 


the momentum relation Eq. (4-5) applied immediately behind the shock wave with pp and yo taken constant and 
equal to their values in the shock layer. Consistently we evaluate the pressure difference p.(R,) — po(R,) from the 
modified shock relations, Eqs. (6-22) and (6-23). The result is 
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We now note that the shock conditions, Eqs. (6-21)- 
(6-25), when taken together with the boundary condi- 
tions at the body, are sufficient to determine a solution 
to Eqs. (5-5)—(5-7). As already mentioned, in arriving 
at the shock boundary conditions, we have neglected a 
conservation equation for h» (or hy) in order to avoid 
overspecifying the boundary conditions. Although the 
viscous correction at the shock involving /2’(R,) is of 
the same order as the other viscous corrections we 
neglect its influence since /2(R,) is by itself small, being 
of the order of h... Furthermore, the function /» is only 
weakly coupled with 4) which is the dominant function 
in the stagnation region where the solution can be 
considered most significant. These remarks are equiv- 
alent to allowing ho! (Rs) to be determined by the 
integration through the viscous layer itself. 

An alternative approach for avoiding overspecifica- 
tion of the modified boundary conditions would have 
been to drop the function h, completely. However, in 
order to match the incipient merged shock boundary 
conditions with the already existing constant-density 
formulation containing h2, the above approximate 
approach was chosen. It must be emphasized again 
that both approximations are equally valid for the 
stagnation region. 

Eqs. (6-21), (6-24), and (6-28), when expressed in 


RiP _ le 4 | Ret,’ | + “| = 
Re,U Re, «eU} U Re, ¢U J 
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terms of V, are the appropriate shock conditions to be 
used for the integration of the vorticity equation, Eq. 
(5-5). One method of carrying out the numerical 
solution of this equation is to choose a value of y and 
Re, and assume a value of «. A value of u’(R,), say, 
is then selected and in the same manner as for the vis- 
cous layer an iteration of this selected quantity is 
carried out until the boundary conditions at the body, 
Eq. (5-8), are both simultaneously satisfied in the 
integration. For any given y, Re,, and ¢, a body loca- 
tion R, can be found. We then proceed to the integra- 
tion of the enthalpy equation for /2, Eq. (5-6). This is 
carried out in the same manner as for the viscous layer 
with the shock condition on fz expressed through Eq. 
(6-27). Eq. (5-7) for ho is then integrated, with not 
only the value of io(R,) but also ho’(R,) initially speci- 
fied, Eqs. (6-22), (6-25), and (6-26). Integration of 
this equation from the shock to the body determines the 
value of ho(R,), which may not agree with the given wall 
value. We then choose a new value of ¢, return to the 
vorticity equation, and repeat the whole procedure 
until the desired wall condition on the function hp is 
obtained. This coupling between the momentum and 
energy equations, despite the continued treatment of 
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Fic. 6. Incipient merged layer surface shear stress on a sphere 
(tem = pPoU? sin 9 cos 8). 


the shock layer as a constant-density region, is to be 
expected because of the introduction of gradients into 
the shock-wave conservation equations. It represents 
essentially a physical coupling between the strength of 
the shock-wave (e) and the body geometry (Re,), de- 
termined by the body temperature (/,), which is not 
present when the shock-wave discontinuity is repre- 
sented by the Hugoniot relations. 

Calculations have been carried out in the manner 
indicated for the incipient merged layer on a highly 
cooled sphere. The same physical constants were em- 
ployed as for the viscous layer calculations, with the 
difference that « couJd no longer be fixed. Here, how- 
ever, the gas is specified equivalently through the isen- 
tropic exponent y which was taken equal to 11/9. For 
Reynolds numbers large enough to be in the viscous 
layer regime, this corresponds through the perfect gas 
strong shock relation Eq. (3-5) to a density ratio ¢ of 
0.1. 

Figs. 6 and 7 show the shear and heat transfer ob- 
tained by these calculations. The most significant 
feature of these results is that, as anticipated, both of 
these quantities drop off relative to the corresponding 
viscous layer solutions, beginning at a Reynolds number 
of about 100. Instead of using Re,, these results are 
plotted using a Reynolds number Re, which is some- 
what more appropriate for the incipient merged layer 
analysis, in which the location of the shock wave is not 
sharply defined and in which values of the fluid proper- 
ties in the viscous layer vary with Reynolds number. 
This Reynolds number Re, is defined by 


Re, = Po UR,/io (6-29) 


where R, is the radius of curvature of the sphere, and 
fi) is the viscosity based on the stagnation enthalpy of 
the free stream, U?/2. All these quantities are speci- 
fied in a flight situation. The relation between Re, 


and Re, is 
Re, = Re,(Ro/Rs) [ho,/ (1/2) U?]*/? (6-30) 


The change from yo, to fp is made by assuming the vis- 
cosity to vary as the square root of the enthalpy on the 
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Fic. 7. Incipient merged layer stagnation-point heat-transfer 
rate on a sphere (—gGpy = (1/2) po U* cos 6) 


stagnation streamline. For the viscous layer calcula- 
tions uo, = Mo because yo, = U?/2. 

Table 2 shows the comparison between the two 
Reynolds numbers, the variation in density ratio across 
the shock wave, and the values of some other important 
flow quantities behind the shock. An important fea- 
ture of these results is the decrease in density ratio and 
enthalpy with decreasing Reynolds number. Note 
that the dimensionless total enthalpy can be obtained 
from the static enthalpy given in the table by adding 
(v,/U)?. 

Calculations have also been carried out in the in- 
cipient merged layer regime for an insulated sphere. 
They were done in the same manner as described above 
for the highly cooled case, except that the enthalpy 
boundary conditions on the body were changed to h,’ = 
ho’ = 0, in accordance with Eq. (4-14). The same 
physical constants were employed. Some results are 
presented in Table 3. In contrast to the highly cooled 
case, the density ratio and enthalpy behind the shock 
increase with decreasing Reynolds number. Note also 
that the stagnation-point enthalpy and the enthalpy 


TABLE 2 
Incipient Merged Layer Results for Highly Cooled Sphere 
y = 11/9, Pr = 0.71, h,/(1/2)U? +0 


y 








ho, Po 

— x a 
Re, € R,/R, (1/2)U? Rep v,/U Pa U? 
100 0.104* 1.10 0.995 90.9 0.097 0.870 
40 0.094 1.10 0.873 34.0 0.076 0.835 
25 0.078 1.08 0.710 19.5 0.048 0.827 


* This differs slightly from the viscous layer value of 0.1. The 
difference may well be attributable to the approximations of the 
theory and small numerical inaccuracies. 





TABLE 3 
Incipient Merged Layer Results for Insulated Sphere 
y = 11/9, Pr = 0.71, mb’ = 0 
he i? & 


Re, e R,/R, (1/2)U2 Re w/U pa? (1/2)U? 
05 36.8 0.097 0.816 1.06 
05 22.8 0.091 0.770 1.09 








40 0.117 1.1 
25 0.124 1.1: 


>— 


t 
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behind the shock are higher than the free-stream stag- 
nation value. Although not shown, the stagnation- 
point pressure behaves similarly. 

We may expect these incipient merged layer results 
to hold down to a Reynolds number of 20 or 30, and 
possibly somewhat lower. To determine this limit 
more precisely would entail calculations based on a 
fully merged layer analysis. 


(7) Discussion of Results and Relation 
to Experiment 


The results of the present analysis clearly show that 
at low Reynolds numbers the usual model of hyper- 
sonic blunt-body flow—an infinitesimally thin shock 
wave, an inviscid region, and a boundary layer—does 
not give a correct description of the flow quantities, at 
least in the stagnation region. The thickening, first of 
the boundary layer and then of the shock wave, as 
Reynolds number decreases, leads to significant changes 
in these quantities. In particular, when the shock 
thickness becomes important, the curvature, viscous, 
and heat conduction effects have a substantial in- 
fluence on the conservation of mass, momentum, and 
energy across the shock wave, as well as in the shock 
layer. This in turn provides a physical coupling be- 
tween shock-wave strength and body size and tempera- 
ture which is not present in the usual flow picture. 

Among the stagnation region phenomena predicted 
by the present theory as the Reynolds number decreases 
are: an increase in detachment distance for insulated 
bodies, an altered behavior of surface shear stress and 
heat-transfer rate relative to extrapolated boundary- 
layer theory, an increase of stagnation enthalpy and 
pressure above their free-stream value on insulated 
bodies, a decrease of normal velocity and enthalpy be- 
hind a shock for highly-cooled bodies, an increase or 
decrease of density behind the shock for highly-cooled 
or insulated bodies, respectively, and a decrease in 
pressure behind the shock wave for highly-cooled bodies. 

Although some of these effects might seem at first 
glance to be unexpected, a number of them have been 
pointed out by other investigators, and in some low 
Reynolds number experiments certain general features 
are found which permit qualitative comparison. 
Moulic” has made afterglow tests on insulated spheres 
at low supersonic speeds. These experiments indicate 
an increased detachment distance with decreasing 
Reynolds number, with the order of the increase con- 
sistent with the present calculations. 

Also of interest are heat-transfer and equilibrium 
temperature measurements on fine cylindrical wires in 
a supersonic stream, for example, those of Laufer and 
McClellan.'* Although the wires were heated and the 
total heat transfer was measured rather than the local 
heat-transfer rate calculated here, the results are none- 


theless significant. The change in heat transfer with 


decreasing Reynolds number, relative to boundary- 
layer theory, is in the direction indicated by our vis- 
caus layer results (Fig. 5). 


Also the measured equilib- 
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rium temperature shows a rapid increase at low Reyn- 
olds numbers, rising above the free-stream stagnation 
temperature. This effect, caused by the added work 
done by viscous stresses, is consistent with the results 
presented in Table 3. 

Then there ate total pressure measurements, for 
example, those of Sherman!® and Enkenhus,” in which 
the pressure is shown to increase with decreasing Reyn- 
olds number at an insulated stagnation point. Again, 
qualitative agreement with the present theory is indi- 
cated. For additional experimental results the reader 
is referred to Schaaf and Chambré.‘ 

While some of our results are therefore found to be 
in qualitative agreement with available experimental 
data, there are a number of other effects mentioned 
above, such as the density, pressure, or enthalpy 
changes, for example, which, if found and measured 
experimentally, would help to indicate whether the 
present analysis correctly predicts the physical phe- 
nomena. (Note added in proof: After the work of 
references 3 and 25, a report by Hoshizaki” appeared 
in which he also considered the constant-density vis- 
cous layer on a sphere. His solutions are not complete 
symmetric ones in the sense of the present paper, and 
the additional approximations made are somewhat more 
restrictive than those used here. However, calculations 
for several values of e were made, with a lower limit of 
approximately 500 for Re,. Where these high Reynolds 
number calculations overlap ours, there is reasonable 
agreement.) 


Appendix (A)—Navier-Stokes Equations in Polar 
Coordinates 


The Navier-Stokes momentum equations, with 
variable viscosity, may be found written in both cylin- 
drical and spherical polar coordinates in reference 
22, pp. 59-61. The two cases may be combined by 
the use of a factor 7, which is zero for the cylindrical 
case and one for the spherical case. The energy equa- 
tion in Cartesian tensor form may be found in the same 
reference, pp. 15-16, and can easily be converted to 
the polar coordinates by the use of the strain tensor 
components given on pp. 59-61. 

We use the notation shown in Fig. 2, with the sub- 
script used to denote the direction normal to the 
R — @ plane in either coordinate system. Since we are 
interested in axisymmetric or two-dimensional flow, 
all derivatives and velocity components in the direction 
of nm vanish. Then the components of the strain tensor 


are 
Crr = 2(Odr ‘OR), 


€re = R(O/OR)(go/R) + (1/R)(OgR/08), 
Cnn = j(2 ‘R) (qe + de cot 6) 


Coo = (2 R)[(Ogo O#) + Ir| 


The components of the shear stress tensor are related 
to those of the strain tensor by the viscosity coefficients 


and are given by 
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Trr = MerR a 2[u* — 3)|¥-q; Teo = Mega + 2[u* - (u/3)|V-q 
Tnn _ Henn + 2[u* ines (u 3)|V-q; TRO and Mere; TRn _ Ten = 0 
where 2V -q = err + C9 + Cnn = (2/R)[R(Ogr/OR) + Gr + (Oge/00) + j(gr + Go cot 6)] 


The momentum equations can now be written as 


op Og Je O”G Je" o 
f = p| ds “ + Pa = Herr + 2(u* — w/3)V-q] + 


oR OR R 00 R OR 
as (were) + . (Cre — Coe) 4 i; (Cre — Cnn t+ Cre Cot 0) 
(uere) + 2 A Cre + J z lere + (€e@ — nn) cot 8] 


The energy equation is 


| oh 4 Jo 4 op a Je Op 4164 l | (er 3) ‘ o (¢ or) a (= 4 cot 6 a) | 
Gr = @D 5 
"Ll oR Rol “aR” Ra RLOR \" aR) * ae\R w/ '”\aR* R ao 
where the dissipation function ® is defined by 
l 2 
m = 9 u(€rr + Cag" + ta + 2€ re") + (u* a =H (va) 


0 


Finally, the continuity equation in these coordinates becomes 


l 0 pc Opge ; 
V-(pq) = al oy + = + pdr + Je(dr + qo cot 0) | = 0) 
Appendix (B)—Shock-Wave Conservation Equations 


The compressible symmetry equations, Eqs. (4-2)—(4-7), are rewritten by eliminating fo’, po’ and pe — po from the 
enthalpy equations, Eqs. (4-6) and (4-7), by means of the momentum equations, Eqs. (4-3)—(4-5). Also, the con- 
tinuity equation, Eq. (4.7), is used for any quantity f to write 


povf’ = (povf)’ — f(pov)’ = (povf)’ — (1 + f)pof(u — v)/R 


With these changes and some rearrangement, Eqs. (4-2)—(4-7) become (with yuo* ignored) 


(pov)’ — (1 + f)po(u — v)/R = 0 (B-1) 
} u—v 4 l+ju- a : Fal “u— "] 
Do’ + (pov?)’ — (1 ) pov =— o{v' +- + (1 +7) u’ — 2v’ — 3 (B-2) 
t (p + j)p R 3 E ( 9 R Dp R 
u—i M8 (w +") ‘a 
Dy — pou = = (B-: 
i p R R R ) 
Mt Real. stint? | ( , &m *] 4 
< = yvu 2 ou = = 3 
R p J)P R Ma R 


I aaa RD yy BED a i 
Pr | (ut ) + R polo + R2 po(he — i) | “+ | molt ot. 2 ) | — 


: 1 \u-v 4 , iw l+ju-—v\]' 
(1 + J) po (i. + 2 e] R — = 3 | mo + 9 R } a 


2 ~~? » Mo : ? ,#—t f x =v ss 
(1 + 7)uv — + (1 + 7) v{u’ — 20° — 3 =~ oe (B-5) 


3 R R 


These equations are integrated on R from R, to R— ~, the boundary conditions at ~, Eqs. (4-15)—(4-17), are in- 
serted, and all gradients are taken to vanish at ~. The result of the integration is then 
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These are the shock-wave conservation equations which are used to obtain the shock-wave boundary conditions in 


the incipient merged layer model of the flow field. 
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A Study of Vortex Cancellation 


ALLAN SCHAFFER* 
Space Technology Laboratories, Inc. 


Summary 


The cancellation of a vortex by means of another concentric 
vortex of equal strength but opposite spin is investigated. When 
such a cancellation occurs, there is a recovery of static pressure. 
The vortices are generated by means of two three-dimensional 
airfoils cantilevered from the duct wall, one being situated in the 
wake of the other. The airfoils have opposite effective angles of 
attack and therefore have trailing vortices of opposite spin, as 
required. 

It is demonstrated experimentally that there exists an optimum 
angle of attack for the second airfoil which cancels the vortex 
from the first airfoil and restores uniform flow downstream of 
the two airfoils. A theoretical solution of this optimum angle of 
attack is presented, and it is found to depend upon the angle of 
attack of the first airfoil and upon the geometrical properties of 
the wings. The pressure recovery accompanying the vortex can- 
cellation is also studied. Theoretical considerations, based on 
the model of a vortex filament in the center of a circular tube 
show that a maximum of 62 per cent of the static pressure drop 
across the first airfoil can be recovered. This maximum is im- 
posed, irrespective of skin friction and separation losses, by the 
irreversibility associated with establishing a vortex field. Ex- 
perimental pressure recoveries of 50 per cent are realized. 

Perhaps the primary value of the present study is the oppor- 
tunity it provides to verify certain of the fundamental concepts 
of fluid mechanics which are brought into play when the trailing 
vortex system of a lifting wing is cancelled by a second wing. 


Symbols 
a = radius of vortex core 
A = aspect ratio 
c = chord 
p = pressure 
Ap = mean duct pressure upstream of airfoil minus mean 
duct pressure downstream of airfoil 
Ap; = irreversible portion of Ap 
Ap, = reversible portion of Ap 
q = velocity head 
r = radial coordinate 
R = duct radius 
$ = half span of the airfoil 
id = half span of the horseshoe vortex 
V = velocity 
a = angle of attack 
B = (A/2) + (2a/r) 
iy = circulation 
U7] = pressure recovery efficiency 
p = mass density 
( ), = upstream airfoil 
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( )e = downstream airfoil 

( )* = quantity evaluated at airfoil midspan 
( )o = undisturbed flow 

( )» = quantity evaluated at the wall 

( )e = tangential component 

(~) = spatial average value 


(1) Introduction 


— PAPER describes a study of the cancellation of 
the trailing vortex system of a lifting wing by the 
use of a second wing situated in the wake of the first. 

Three-dimensional half-wings cantilevered from the 
duct wall are used. A single rolled-up vortex emanates 
from each. The vortex wake downstream from a half- 
wing consists of a relatively thin filament, in which to 
a good approximation the curl of the velocity is con- 
stant, surrounded by a flow in which the curl of the 
velocity is zero. It follows that the tangential veloc- 
ity varies directly with radius within the filament and 
inversely with radius outside the filament. Hence, the 
tangential kinetic energy of the vortex flowfield emanat- 
ing from the upstream wing is in a highly ordered form 
and is thus subject to recovery. If the flow can be 
straightened and the vortex cancelled without appre- 
ciable increase in entropy, the tangential kinetic energy 
can be recovered in the form of increased static pressure. 

Pressure recovery by this technique is of interest in 
the design of rotating compressors; more novel ap- 
plications are also conceivable. For example, it is 
possible to mix fluids by means of a vortex, then cancel 
the vortex and recover a substantial portion of the 
energy assigned to do the mixing. The net result is a 
mixing chamber with low overall pressure drop. How- 
ever, it is not the purpose of this paper to discuss prac- 
tical applications of vortex cancellation but rather to 
present the results of some fundamental experiments 
and theoretical analyses concerning the phenomenon. 

Experimental studies of vortex cancellation have been 
carried out in a conventional low-speed wind tunnel. 
The technique involves first the generation of a vortex 
by a half-wing with a given angle of attack; then ata 
point downstream a similar wing with an opposite 
effective angle of attack generates a vortex of equal 
strength but opposite sense of rotation. The wings are 
aligned so that the two opposite-spinning vortices are 
concentric, and in this manner vortex cancellation can 
be achieved. 

For maximum cancellation, the angle of attack of the 
second airfoil must be optimized so that the strength 
of its trailing vortex is just equal in magnitude to that 
of the first airfoil. This optimization problem is inves- 
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Fic. 1. 6 by 6 in. test section. 


tigated experimentally and theoretically in the paper. 
Similar studies are made of the pressure recovery ac- 
complished by the second airfoil. 


(2) Description of Experimental Facilities 


The wind tunnel is of conventional low-speed design. 
It has a maximum test section velocity of 200 ft. per 
sec. and a 6 by 6 in. square test section, as shown in 
Figs. 1 and 2. The wings are cantilevered from 
pivotable discs which form part of the tunnel walls 
and which can be rotated to any convenient angle of 
attack. The centers of the two airfoils are 11 in. apart. 
A tuft grid, located 5 in. past the center of the down- 
stream airfoil, is used to observe the vortices.' Static 
pressure measurements are made by means of both wall 
taps and spherically ended probes. 


(3) Experimental Techniques 


The two airfoils are usually geometrically identical. 
The angle of attack a; of the upstream airfoil is gener- 
ally set just under the stall point so that the strongest 
possible vortex without separation is produced. All of 
the airfoils tested are of relatively low aspect ratio (in 
the range 1 to 4) so that the trailing vortex sheet rolls 
up into a single vortex in a distance downstream of the 
order of one chord’; for this reason no airfoil pairs 
have been placed closer together than two chords in 
the present experimental study. 

In order to achieve full vortex cancellation, the fila- 
ments from the two airfoils must be coincident. Since 
the filaments exhibit a neutral stability in that they 
neither attract nor repel each other, rather careful posi- 
tioning of the wings is necessary to ensure that the 
filaments are precisely concentric. As predicted by 
vortex theory, a slight mismatch still results in cancella- 
tion of the irrotational field, but the rotational cores do 
not then cancel; instead, they translate each other 
violently. 

Once the airfoils have been aligned so that the vor- 
tices are concentric, the angle of attack a2 of the down- 
stream airfoil is optimized by rotating the airfoil mount- 
ing disc until the angle is found at which the minimum 
ordered rotation is visible in the tuft grid pattern. 
The uncertainty in determining the optimum az by this 
technique is about +1°. 


(4) Photographic Results of Vortex Cancellation 
Experiments 


Cancellation tests have been made for three sets of 
airfoil pairs, the specifications for which are given in 
Table 1. All are NACA 0012 symmetrical untwisted 
half-wings of rectangular plan form. Test results are 
reported for Vy) = 196 ft. per sec.; however the phenom- 
enon exhibits negligible Reynolds number dependence 
within the range of the tests. 

Figs. 3 and 4 are typical tuft-grid photographs of the 
cancellation phenomenon for airfoil pairs (2) and (3) 
of Table 1, respectively. In Figs. 3(a) and 4(a) only 
the upstream airfoil is installed, whereas in Figs. 3(b) 
and 4(b) the downstream airfoil is also installed and 
ais optimized. Depending upon the illumination, the 
airfoils are sometimes visible in the photographs. As 
seen in the figures, the cancellation is virtually com- 
plete except for some slight residual core vorticity. 
The optimized values for az are —3!/2° + 1° with air- 
foil pair (2) and —2!'/.° + 1° with airfoil pair (3). 

Although only the cancellation phenomenon is shown 
in the pictures, it is simple to obtain other interactions 
such as vortex addition merely by varying az. 


(5) Theoretical Analysis of Vortex Cancellation 


The present analysis is intended to be only a first- 
order description of the phenomenon. The model in 
Fig. 5 is assumed. The upstream airfoil is represented 
by a horseshoe vortex with strength Ty. The bound 
portion of the horseshoe vortex has negligible effect 
and is ignored in the derivation. The discussion is 
limited to symmetrical untwisted wings with the up- 
stream and downstream airfoils being geometrically 
identical. 

The basis of this analysis is that the optimum ay 
will occur when a trailing vortex from the downstream 
airfoil has strength of equal magnitude and opposite 
sign to the corresponding vortex from the upstream 
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A STUDY OF 


VORTEX CANCELLATION 








TABLE 1 


Comparison Between Experimental and Theoretical Cancellation Results 


Pair No. r id A 
(1) 13/, in. 31/4 in. 4.00 
(2) 3 in. 23/, in. 1.94 
(3) 5!/ in. 31/3 in. 1.16 


airfoil. In view of the fact that the circulation Ty 
at the midspan of an airfoil equals the strength of one 
of its fully rolled-up vortices, the requirement can be 
formulated as Tx. = —Ix:. It is logical that, for this 
requirement to be met, a, should equal the negative 
of a, decreased in magnitude by the induced angle of 
attack resulting from the upwashf generated by the up- 


stream airfoil. In equation form, one obtains 


a, = —a, + (I'e1/7s'Vo) (5.1) 
Eq. (5.1) can also be derived on a more formal basis 
by means of elliptical wing theory. 

In proceeding with the analysis, the theory for ellip- 
tically loaded wings of small span, developed by 
Weinig,* will be used, since the wings employed experi- 
mentally have relatively low aspect ratios. The re- 
sults given by conventional lifting-line theory can 
always be obtained as a special case by examining the 
asymptotic limit of the formulas for large aspect ratio. 

A useful equation yielded by Weinig’s theory is 


Tx = mex Vo [8 tan B-'/(1 + tan B~") Ja (5.2) 


where 6 = (4/2) + (2a/r). Combining Eq. (5.1) 
with Eq. (5.2) evaluated for the upstream airfoil, one 
obtains the desired relationship between a2 and ai; 
Viz., 


a. = } [cx tan B-!/s’ (1 + tanB-!)] — lfa (5.3) 


As might be expected, no velocity dependence is ob- 
tained. For A > 4a,;/z, the formula reduces to that 
given by conventional lifting-line theory; viz., 


a. = | [cx/s’(1 + 2/A)| — lb (5.4) 


Eq. (5.4) shows that for wings with higher aspect ratio, 
a varies directly with a. 

In Table 1 the predictions of Eq. (5.3) are compared 
to the three experimental data points available, all of 
which are for rectangular wings. The experimental 
uncertainty in determining a: is about +1°. There 
appears to be a first-order agreement between theory 
atid experiment; however, no final conclusion as to 
the accuracy of the theory can be made because insuffi- 


cient experimental data are available. 


(6) Pressure Recovery by Vortex Cancellation 


The pressure drop Ap across an airfoil will be de- 
fined as the mean duct pressure upstream of the wing 
minus the mean duct pressure downstream. More- 


+ The term ‘‘upwash”’ is used because the interaction is such 
as to increase the angle of attack of the downstream airfoil. 


stall 


13° 








Q2)exp. —~ G2 /theory 
a Q2 Jexy Q2 )theory a 
10° — 5° —6.2° 0.12 
14° —3'/,° —4.5° 0.07 
18° —2!/,° 1.0° —0.19 


over, Ap can be divided into a reversible portion Ap, 
associated with establishing the ordered vortex field, 
and an irreversible portion Ap;, associated with mech- 
anisms producing disordered motion. When the vortex 


is cancelled, Ap,» = —Apy. In view of this relation- 
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Fic. 3. Cancellation for airfoil pair (2), Vo = 196 ft. per sec. 
(a) Left: one airfoil a = +14°. (b) Right: two airfoils, a, = 
+14°, a= —3! gs! 
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Fic. 4. Cancellation for airfoil pair (3), Vo = 196 ft. per sec. 
(a) Left: one airfoil, a = +18°. (b) Right: two airfoils, 
a, = +18°, a2 = —2!/,° 
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Fic. 6. Tangential velocity distribution. 
V, = 1T,/2ra?, ra 
Vo = Ty /2zr7, esrsR 


ship, a recovery efficiency 7 can be defined as 

n= 1-— [(Apa + Api)/Api] (6.1) 
To simplify the computation, the assumption is made 
that 


Apa = Api (6.2) 


This assumption is reasonable since the airfoils are 
geometrically identical and since for cancellation they 
experience much the same relative flows. Eq. (6.2), 
when substituted into Eq. (6.1), together with the 
identity Ap; = Apn + Apu, yields 


2=1- {2 {1 + (Apa /Apa) }} 


Eq. (6.3) is useful in determining both the experi- 
mental efficiency and the maximum theoretical effi- 
ciency. Experimental efficiencies determined this way 
obviate the tedious task of making pressure traverses 
in the wake of an airfoil. The procedure is first to 
ascertain the vortex strength Ty and core radius a in 
the wake from the first airfoil by means of a tuft grid 
photograph. Assuming the usual velocity distribution 
for a vortex, Fig. 6, one can computef Ap, knowing 
Ty and a. Then Ap, is found by measuring the total 
static pressure drop across the two airfoils with vortex 
cancellation.{ By Eq. (6.2) Apa equals 1/2 this value. 
Thus, knowing Af, and Api, one computes the ex- 
perimental » from Eq. (6.3). For the airfoils tested, 
the value obtained by this method is approximately 
n = 0.5, and a similar value results when 7 is determined 
by means of pressure traverses. Thus, experimental 
pressure recoveries of approximately 50 per cent have 
been realized. 


(6.3) 


t See Eq. (6.6). 
t Wall taps can be used for this measurement since there is 
uniform flow upstream of the first airfoil and downstream of the 


second one. 
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To determine the maximum theoretical value of », 
consider steady incompressible flow over a half airfoil 
in a circular duct of constant area with no wall friction. 
We must compute Ap, and Ap; in Eq. (6.3). Subscript 
( )o denotes a station in the undisturbed flow upstream 
of the first airfoil, no subscript denotes a station in 
the wake shown in Fig. 6, and subscript ( ), denotes the 
tangential component in the wake. From an energy 
balance together with the relationship V? = V,? + 
V,” it can easily be shown that, for the upstream air- 
foil, 


Ap; = Po = p “i do (6.4) 


where the notation (~) means average value over the 
cross section and g, = (1/2) p V,? is the tangential 
Ap, is by definition the value of Ap 
Setting Ap; = 0 in Eq. (6.4), one 


velocity head. 
when Ap; is zero. 
obtains 


(6.5) 


bo — Pp = % = AP, 
Eq. (6.5) can be evaluated from the velocity profile 
equations—namely, V, = 7I'x/2ma? forr < aand V, = 
T'y/2mr fora <r<R. One obtains the final expression 
for Ap, as 


Ap, = 2 dow[In(R/a) + (1/4)] 


where Qo» is the tangential velocity head at the wall. 

It remains to compute Ap;. The minimum value 
of Ap; can be determined by specifying the skin fric- 
tion and separation losses to be zero. The residual 
Api, called Ap;)min, results from the irreversibility asso- 
ciated with establishing the vortex and is evaluated by 
means of Eq. (6.4). The value of g, to be used in Eq. 
(6.4) is given by Eq. (6.5) together with Eq. (6.6). 
The expression for / can be deduced by applying New- 
ton’s Second Law along the radial direction in Fig. 6; 
this gives 


P = Pw — 2 dow[In (R/a) + (1/4)] 


(6.6) 


(6.7) 
Combining Eqs. (6.5), (6.6), and (6.7) with Eq. (6.4), 


one obtains 
AP) min = Po = Pw (6.8) 


But Bernoulli’s equation is applicable between points 
in the undisturbed flow and points in the irrotational 
region downstream of the airfoil, so that 


(Po ‘p) + (Vo?/2) = (DPw/p) + (Vy2/2) 


and since V,? = Vo? + Vow”, Eq. (6.8) reduces to 


(6.9) 


APi) min = Vow (6.10) 


When Egs. (6.6) and (6.10) are substituted into Eq. 
(6.3), the theoretical maximum value of 7 is obtained 
as 


Nmazr = [41n (R/a) — 1]/[41n (R/a) + 3] (6.11) 


It is significant that ymez is a function only of R/a and 

not even a very strong function of it. The numerical 

resuits differ little whether R/a or s/a is used in Eq. 
( Continued on page 233) 
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Optimization of the Flight-Control, Airframe 
System! 


I. L. ASAHKENAS* ano D. T. MCRUER** 


Systems Technology, Inc. 


Summary 


The analysis of a system composed of a flight control and an 
airframe is ordinarily using servoanalysis 
techniques in which the airframe and controller are described in 
Optimization of such a system at 


accomplished by 


terms of transfer functions. 
the analytical level involves the selection of desirable relative 
locations for the poles and zeros of the airframe and controller 
transfer functions. The appropriate airframe pole and zero 
locations thus selected are then reflected back into the airframe 
geometrical characteristics. To achieve such an optimum sys- 
tem thus demands an overall approach which considers both the 
airframe and flight controller to be alterable system elements. 
This paper develops the optimization concept and illustrates the 
procedure with practical examples. 


Symbols 
a,A = polynomial coefficients 
ay = lateral acceleration at the center of gravity 
ie = lateral acceleration at a distance /, from the c.g., 
ay’ = ay + let 
a, = normal acceleration at the center of gravity 
a,’ = normal acceleration at a distance /, from the c.g., 


; : 
a,’ = az — 12g 
AR = amplitude ratio 


A,(s) = transfer function relating lateral acceleration to 
control-surface deflection, a,(s)/5(s) 

A,,‘(s) = transfer function relating lateral acceleration at a 
distance /, from the c.g. to control-surface de- 
flection, a,’(s)/5(s) 

1,(s) = transfer function relating normal acceleration to 
control-surface deflection, a,(s)/5(s) 

b = wing span 

b, B = polynomial coefficients 

et = polynomial coefficients 

C1 = rolling-moment coefficient, 

(Roll Moment)/(1/2)p Uo?.Sb 

Ci, = dihedral parameter, variation of rolling-moment 
coefficient with sideslip, 0C;/08 

Cis = aileron rolling effectiveness, 0C:/05¢ 

a 

Cy = roll-damping coefficient, 0C;/0(pb/2 Uo) 

Ce = roll coefficient due to yawing velocity, 
0C1/2(rb/2 Uo) 

Cr = lift coefficient, nW/(1/2)pUo2S 

CG, = yawing-moment coefficient, 


(Yawing Moment)/(1/2)pUo2.Sb 


Presented at the Aerodynamics—III Session, IAS 27th Annual 
Meeting, N.Y., Jan. 26-29, 1959. Revised and received July 
20, 1959. 

+ This research was supported in whole or in part by the 
United States Air Force under Contract No. AF 33(616)-3906, 
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Wright Air Development Center, Wright-Patterson Air Force 
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Ms 
M, 
M, 
Me 
Me 
M.' 


Na,(s) 
Na y “(s) 


Ng 
Nas) 


Ns 


static directional stability, 0C,/08 
aileron yawing effectiveness, O0C,/0d¢ 


OC,,/2( pb/2 Uo) 
OC, /0(rb/2 Us) 
lateral-force coefficient, 
(Lateral Force)/(1/2)pU,2S 
variation of C, with sideslip angle, 0C,/0p 
polynomial coefficients 
abbreviation for decibels = 20 logy | 
polynomial coefficients 
error 
acceleration due to gravity 
altitude 
transfer function relating 
surface deflection, h(s)/5(s) 
moments of inertia referred to stability axes 
product of inertia referred to stability axes 
the imaginary portion of the complex variable 


altitude to control- 


s=o2tjw 

radii of gyration about X, Y, and Z axis, respec- 
tively—e.g., Jz = mk,? 

transfer function of ideal simplified yaw-rate con- 
troller, neglecting servo and sensor dynamics 
and assuming a linear system, 6,/r 

transfer function of ideal, simplified roll controller 

the zero-frequency value of the transfer function 
5, 

distance along the X axis from the c.g., positive 
forward 

pSUbCi, Rie 

ft 

pSU9*bCi,/21z 

pSL 0b? Ci,/41z 

pS Uob?Ci,/41z 

Lg/ Uo; pS bCig | 

[Li + (122/Iz)Ni) 


to any motion, or input, quantity 


{1 — (/,,?/I,J,)]| where 7 refers 


Laplace transform of 

mass 

UM. 

UMu 

pSUv*cCu, 2I, 

pSl ne 2Cay ‘41, 

pSL f cCm,/ Ty 

pSUgeCu,/21y 

pSc?Cu,,/41y 

M; — 2Z;:M;/Z;, where 7 refers to any motion, or 
input, quantity 

numerator of transfer function 
lateral acceleration to control deflection 

numerator of transfer function relating lateral 
acceleration (at point /,) to control deflection 

pSUo*bCng 2I, 

numerator of transfer function relating side lip 


relating c.g. 


angle to control deflection 
pSU,*bC,;/21, 
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Nya(s) 


Nias) 


Nels) 


Nos) 
N(s) 


Nis) 


Ras) 


L/T; 


Uo 
Uxs) 


ll 


Il 
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numerator of transfer function relating directly 
controlled airframe-motion quantity to dis- 
turbance input 
numerator of transfer function relating directly 
controlled airframe-motion quantity to control 
deflection 
numerator of transfer function relating indirectly 
controlled airframe-motion quantity to control 
deflection 
numerator of transfer function relating indirectly 
controlled airframe-motion quantity to disturb- 
ance input 
numerator of transfer function relating pitch angle 
to control deflection 
numerator of transfer function relating roll angle 
to control deflection 
numerator of transfer function relating rate of 
climb to control deflection 
pS Uob?Cn,/41s 
pS Uob?C,,,/41- 
numerator of transfer function relating yaw rate 
to control deflection 
numerator of transfer function relating airspeed to 
control deflection 
Nep/ Uo 
numerator of transfer function relating vertical 
velocity to control deflection 
numerator of transfer function relating lateral 
angular velocity to control deflection 
[Ni + (Le2/Ie)Lil/(1 — (I22?/I2I2)], where i refers 
to any motion, or input, quantity 
roll rate, angular velocity about the X axis, posi- 
tive right wing going down 
pitch rate, angular velocity about the Y axis, posi- 
tive nose going up 
yaw rate, angular velocity about the Z axis, posi- 
tive nose going right 
transfer function relating yaw rate to control- 
surface deflection r(s)/5(s) 
Laplace operator, o + jw 
inverse time constant, particularized by the sub- 
script 
linear, perturbed velocity along the X axis 
linear, steady-state velocity along the X axis 
transfer function relating airspeed to control de- 
flection, u(s)/6(s) 
linear, perturbed velocity along the Y axis 
linear, perturbed velocity along the Z axis 
transfer function relating vertical velocity to con- 
trol deflection, w(s)/5(s) 
UpXw; pSU0%( Cr — Cr,)/2m 
pS Uo —Cp —- Cp,)/m 
pSU( Ci — Co,)/2m 
transfer function of controller forward-loop ampli- 
fication and equalization 
transfer function of sensor- and feedback-loop 
equalization 
open-loop transfer function; the product of all 
series transfer functions in the control loop 
transfer functions with s = s, s = jw, and s = 
—o, respectively 
pSUpCy,/2m 
pSU;Cy,/4m 
pSUoC,,/4m 
pSUCy,/2m 
Zw 
pSUo% —Cx;)/2m 
pS Uoe( — Crq)/4m 
pSU(—-CL — Cr,,)/m 
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aa pSU(—Cr, — Cp)/2m 
Zw = pSce(—Cr,)/4m 
a = instantaneous angle of attack, under no wind 
conditions, w/ U» 
B sideslip angle, under no wind conditions, v/ Uo 
Ba(s) = transfer function relating sideslip to control sur- 
face deflection, B(s)/5(s) 
v = perturbed flight-path angle, 6 — a 
Yo = flight-path angle, steady state 
6 = control-surface deflection 
A = increment change 
A(s) = denominator of airframe transfer functions; char 
acteristic equation when set equal to zero 
Atat. = denominator of the lateral transfer functions; 
lateral characteristic equation when set equal 
to zero 
Atong- = denominator of the longitudinal transfer functions; 
longitudinal chara-teristic equation when set 
equal to zero 
oi = damping ratio of linear second-order system 
particularized by the subscript 
n = radius of the circle generated by the real poles and 
zero (—1/Ty, —1/Tr, —1/T,) on the root- 
locus plot 
0 = pitch angle; fg dt in straight and level flight 
Or = gyro tilt angle measured from the XY axis in the 
X-—Z plane, positive nose up 
04(s) = transfer function relating pitch to control-surface 
deflection, 6(s)/5(s) 
ny = the total directly controlled airframe-motion 
quantity 
Ag = transfer function relating general directly con- 
trolled airframe-motion quantity to control 
deflection, \(s)/5(s) 
Be general disturbance input 
v = distance, on the root locus plot, from the zero, 
—1/T,, to the nearest complex zero or pole 
3 = the total indirectly controlled airframe-motion 
quantity 
8 = transfer function relating indirectly controlled 
airframe-motion quantity to control deflection, 
&(s)/6(s) 
p mass density of air 
= the real portion of the complex variable s = ¢ + jw 
roll angle; (cos yo Pi pdt — sin yo Di rdt) in straight 
and level flight 
4(s) = transfer function relating roll to control surface 
deflection, g(s)/8(s) 
Wi = undamped natural frequency of a second-order 
mode, particularized by the subscript 
w = lateral angular velocity sensed by a gyro tilted 
amount 67, # = rcos@7 + psin 67 
3(s) = transfer function relating lateral angular velocity 
to control-surface deflection, w(s)/6(s) 
Special Subscripts 
a = aileron 
d = Dutch roll 
l = lateral long period 
p = phugoid 
r = rudder 
R = roll subsidence 
s = spiral 
Ss = short period 
OL = open loop 
CL closed loop 


Introduction 


ib gareranges in an engineering sense, involves a 
selection between alternatives within prescribed 
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OPTIMIZATION OF THE 
rules. Under a given set of circumstances, the rules 
or optimizing criteria may be subject to individual dis- 
cretion, judgment, or instinct. The optimization 
process, regardless of the particular rules employed, will 
be far less subjective since it involves the establishment 
of several competing systems based on physical reality. 
The quality of these competing solutions, and the ulti- 
mate success of the optimizing effort will depend di- 
rectly on the degree of understanding which can be 
brought to bear on the various interactions involved 
in the physical system. 

The objectives of the present paper are to lay some 
groundwork for a better appreciation of the overall 
problem of integrating the flight-control, airframe 
system to allow consideration of the airframe as an 
alterable element in the system during preliminary de- 
sign, to discover what is important in the overall sys- 
tem, and finally to determine what can be done to 
simplify the system. The possibility of realizing sub- 
stantial benefits in reduced weight or cost, or in in- 
creased reliability, is the ultimate goal. The approach 
makes use of servo techniques because of their power 
and applicability, but attention is mainly focused on 
the airframe and its reduction to simplified approxi- 
mate servo terms. In outline, the paper starts with a 
presentation of the airframe equations of motion and 
progresses rapidly through the reduction of these 
equations to approximate, factored transfer functions. 


Longitudinal Set: 
(s— X,)u — (Xzs + Xy)wt+ 


—Z,u + (s — Zs — Zy)w + 


—-Mu —-(M,zs+M,)w + 
Lateral Set: 
(s— Y,)v— et == $+ 
‘ 
—L,v + (s — L,’)p 
—N,v — N,'p 


In these equations, the primed lateral derivatives are 
such as to eliminate the explicit appearance of product- 
of-inertia coupling terms between roll and yaw (see, for 
example, references 2 and 3); and the control-surface 
deflection, 5, for the longitudinal set is considered to 
be an elevator and no subscript is used. For the 
lateral set, subscripts a and 7 are later used to dis- 
tinguish between aileron and rudder. 

All the above equations have been written in terms of 
fundamental airframe motion quantities (angular and 
translational velocities) measured about and along an 
axis system fixed in the aircraft; the angles 6 and 9, 
which orient the gravity vector to the airframe axis 
system, have been replaced by their equivalents for 
straight and level flight—i.e., 
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These define the airframe in servo terms and allow a 
ready synthesis of the complete system. A simplified 
example is then used to illustrate possible effects of air- 
frame alterations upon the complete airframe-controller 
system. These results are summarized and _ inter- 
preted to support the major thesis of this paper that an 
overall approach which considers the airframe, as well as 
the control system, to be alterable helps substantially in 
the establishment of the competing systems essential 


to optimization. 


Airframe Transfer Functions 


When the airframe is considered as an element in an 
automatic control system, its properties are defined by 
a series of transfer functions which relate transformed 
output quantities (various airframe motions) to input 
variables (control deflections or external disturbances). 
The transfer functions are ratios of higher-order poly- 
nominals in the Laplace transform variable, s, with 
coefficients composed of stability derivatives and in- 
ertial and gravitational parameters. These forms are 
all well known (see, for example, reference 1), so their 
development is simply summarized below. 

The conventional linearized rigid-body equations 
which describe the perturbed lateral and longitudinal 
airframe motions about straight and level operating 
conditions are given in Laplace transform style in Eqs. 
(1) and (2); initial conditions are assumed zero, and 
inputs due only to control deflections are considered. 


—X,s + g cos Yo 


q = X 56 
Ss 
(Uo + Z,)s + g Sin Yo i (1) 
q = Z 50 
5 
(s— M,)q = M,6 
s - 
—L,'r = L,'8 | (2) 
+(s— N,’)r = N,'5 
@ = £(fqdt) = q/s j 
¢ = L£(cos yo J pdt — sin yo Sirdty (3) 


= cos yo(p/s) — sin yo(r/s) 


The conventional substitutions for angle of attack 
(a = w/U>) and sideslip (8 = v/Uo) can also be made 
in Eqs. (1) and (2), since the assumption of zero air- 
mass velocity is implicit in these equations. 

In addition to the fundamental motion quantities 
given above, many flight control problems also require 
altitude, normal and lateral acceleration (both at the 
center of gravity and at a location /, feet in front of 
the c.g.), and the lateral angular velocity in the X-Z 





200 JOURNAL OF THE AERO/SPACE SCIENCES—MARCH, 


plane displaced through an angle 67 from the X axis. 
These quantities are given below. 


h = Uovo + You + Uny 
a, = w— Ug+t+ gsin yo JS qdt 
= Z,hw+ Zu +26 + 249 + Ziv 
a,’ = a, — 1.4 | 
a, = 0+ User — gsin yo frdt — } (4) 
gcos yo J pdt 
= V5.6, + Vsada + Yop + Vir + Yr 
ay’ = a,t+1e 
w = 7cos@r+ psin Oy J 


The values of some derivatives are such that their 
contributions to transfer function coefficients are 
usually negligible. These include X;, Xi, X4, Zi, Z,, 
Y,*, and Y,*. All of these are assumed to be zero. 
Finally it should be observed that the detailed dis- 
cussions in this paper are restricted to horizontal flight 
(vo = 0) which is the simplest of the three degrees of 
freedom straight and level conditions. With these 
restrictions, the transfer functions of interest become 


Longitudinal Set: 


O(s) _ Nels) _ Ast + Bes + Ce 


1960 


w(s) a No(s) “7 A,s* + Bys*? + Cos + D.. 


W;(s) = - = 
. } 6(s) Dioag. Aiong. 
(6) 
; u(s) N,(s) A,s? + Bys + C, “ 
U;(s) = - = = - (4) 
6(s) Bienes. A ong. 
h(s) N,(s) Ays* + Bys? + Cis + Dp 
Has) = —— = - - = - = 
6(s) SAitsie. SBiteux. 
(8) 
A.3(s) = a.(s)/6(s) = —[sN;,(s)/Atonz.] (9) 
where Ajong. = AS4‘ + Bs? + Cs? + Ds +E (10) 
and 
A= j 
B= —(M, + Xx + Zw + M3) 
C= MZ, -—-Mi,+ X.(M, + Ze + Mz) - 
A+ (2%) 
D= -X,(M.Z, — M,) - 


MiXat MXuZu + e(M3Z, + M,) 
E = g(MyZy — M,Zw) 


The numerator coefficients are summarized in Table 1. 
It should be noted that derivatives with respect to 
both a and w, as well as primed longitudinal derivatives, 



































Q;(s) = -—~- = (5) 
6\° e 2 e e ° 
6(s) Citcus. Pitone:- defined in the Symbols, are used to achieve conciseness. 
TABLE 1 
Summary of Longitudinal Transfer Function Numerators 
A B Cc D 
Exact 2M, + My 2, (M,, - MX) 2, (MX - MX) 
“ee +z) \ 
8 * be ke . My (Z.X,, A X7.,/ 
Approx. M, Z.M, - MaZ 
w 2, - 2 (M, + x.) + UM, X (25Mq - U,Ms) (2M - MZ. 
Exact - 2, (mM, +MX ) 
Ww qw eZ M! 
+ M.(X, - 6) e(Z ‘cn ) os oa 
u 25x Me, es > ws wd U 
Approx. Ms (X, - g) 
Exact 
25 [ My - XM, + m,)| - MZ 2, (MX, -@™M -MX 
~Ri(izext « a 
h - 2 Zs (Ms + M, - x.) ua u x 
Approx. ' > ut (y -\ut 
2M) - 2 [x 7 (xX, - g) “| 
- Zy (My sia “Ud x,) if ES - SA, X, (My , ",)] 2, (MX, - MX, + AM) 
c = I t 7 - - h - @ M 7 ¥ -~ + ~~ (aa 
Exact Z, 1 (2M, + My 1, [ 254, MX.) ay [20M MX) + M, (ZX, ZX Z.) 
“2 My (Z,, si x)| i My (2%, ‘ x2) | 
a 
Zz 
Zz -1) - Z.(Me #M_ + X_) 
5 
Approx. x te i Malis - 2M! > 12.(XM! - XM‘) Z| XM! - (x_ - g)M! 
i Re (2 +x) a xo. wu uw te) ua x u 
x 25 Ww Ms Ww u 
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Lateral Set: 


B(s) N,(s) A,gs* + Bags? + Cys + Dg 


B(s) = = > 
( 6(s) Aiat. Aint. 
(12) 
: N,(s) As? Bs é.. 
®,(s) wal ¢(s) = lS = ¢ = ¢> i ¢ (13) 
6(s) Anat. Anat. 
r(s) N,(s) A,s* + Bs? + Cs + D, 
R;(s) = = = 
6(s) Anat. Atat. 
(14) 
a,(s) Nu, (S) 
A sSj= 2 - = y = 
wae oar 
Aays* + Bays' + Cus? + Days + Ey (5) 
2) 
Atat. 
: , Ss N, d AY 
Ays'(s) oy (s) = v (s) = 
6(s) Atat. 
Aay’s* + Buy'st + Cu/s? + Dy’ + Bay oy 
Atat. 
w(s) N,,(s) 
Q(s) = = - 
als) 5(s) Atat. 
(A,s? + Bs? + C,s + D,) cos 67 (17) 
4 
Aiat. 
where Atat. = ast + bs? + cs? + ds +e (18) 


When the primed derivatives are used, coefficients of 
Ajat. are given by 


a= 1 
b —(Y, + L,’ + N,’) 
c N,’ + L,’(Y, + N,’) — N,’L,’ + Y,N,’ 
d = N,’L,’ — L,'N,’ + 

Y,(N,’L,’ — L,’N,’) — (gLg’/ Uo) 
e = (g/Uo)(L,/N,’ — N,’L,’) 


I 


(19) 


The numerator coefficients are summarized in Table 2. 
Unprimed derivatives are used in this table since they 
lead to slightly simpler results for the numerator factors. 
These can be changed to primed derivative form by the 
simple process of setting 7,, equal to zero and replacing 
the L;, and N; by L,’ and N,’. It should be noted that 
whereas the factors of a given polynomial will be inde- 
pendent of the use of the prime or unprimed derivatives, 
the transfer-function gain must be computed using the 
same type in both numerator and denominator. 


Airframe Transfer Functions in Approximate 
Factored Form 


For maximum utility it is desirable to have the nu- 
merator and denominator polynomials in factored form. 
Each transfer function will then be made up of a ratio 
of first- and second-order polynomial products, and a 
gain constant. The roots of the denominator are the 
transfer function ‘‘poles’’; those of the numerator are 
its “‘zeros.’”’ The gains, poles, and zeros define the 


fundamental properties of dynamic elements and are 
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essential in most servoanalysis and synthesis methods 
and in response calculations. The factorization re- 
quired to find the poles and zeros can be time-consuming 
since most of the polynomials involved are of third or 
higher order. More important, however, is the fact 
that the poles and zeros of these higher order poly- 
nomials are difficult to define in terms of literal values 
of the transfer function coefficients or of the stability 
parameters. This is unfortunate because a large 
number of numerical calculations are then required for 
comprehension of the tie between configuration char- 
acteristics (the stability derivatives, etc.), and the air- 
frame transfer and response characteristics. 

To alleviate this situation, it is highly desirable to 
establish approximate expressions for the poles and 
zeros in terms of the stability derivatives and inertial 
parameters. This approach has often been used in the 
past, either directly by decomposing a fairly complete 
transfer function form, or indirectly by the use of fewer 
airframe degrees of freedom. Both schemes have been 
moderately successful in some instances, while in others 
the results have been severely restricted in application 
and scope. In most cases approximate factors have 
been developed in an artistic fashion, with ill-defined 
conditions under which they are valid. Such pitfalls 
have been skirted as much as practically possible in the 
derivation of the approximate factors given in this 
paper. Where the basic polynomials are of second 
or third order, the approximate factors are derived by 
the use of a method which can be made exact. The 
method is based upon servoanalysis of an “‘equivalent”’ 
closed-loop system, and the degree of approximation in 
any given case can be readily assessed. The approxi- 
mations for the roots of either the lateral or longitudinal 
characteristic equations, which are both fourth order, 
are not established in such a formal manner. Instead, 
they are based upon simplifying assumptions regarding 
the relative magnitudes of various combinations of 
stability derivatives. The validity of these simplifying 
assumptions has been demonstrated, in that the result- 
ing approximations have been applied successfully to a 
large number of actual configurations. 

Table 3 summarizes all the approximate factors for 
the longitudinal case. Table 4 presents the lateral de- 
nominator factors, Table 5 presents the numerator terms 
for aileron inputs, and Table 6 presents the rudder 
numerators. Conditions of validity are shown for all 
cases. The details underlying the derivations pre- 
sented in these summary tables are given in reference 3. 

The ready availability of the relationships given in 
Tables 3-6 is essential to the optimization thesis ad- 
vanced here. They are especially helpful in (1) de- 
veloping the insight required for the determination of 
airframe-autopilot combinations which offer possible 
reductions in overall system complexity, (2) assessing 
the effects of configuration changes on aircraft response 
and on airframe-autopilot system response, and (3) 
showing the detailed effects of particular stability 
derivatives (and their estimated accuracies) upon the 
poles and zeros, and hence upon aircraft and airframe- 
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autopilot responses. 
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The factors are 


OF THE 


also useful in 


dealing with problems involving the airframe alone, such 
as obtaining stability derivatives from flight test data. 


The Controller as an Equivalent-Airframe 
Element 


The foregoing discussion and the approximations 


presented in the tables show that a close tie exists be- 


tween given airframe stability derivatives and their net 


dynamic effect. 


From this point of departure it is 


natural to examine the influence of a control loop by 


considering its role as an ‘‘augmentor’’ of existing 


stability derivatives or a 


sé 


creator” of new ones. 


This 


technique, which is in airframe terms rather than in 
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those of servoanalysis, is called the ‘‘equivalent-deriv 
It is the simplest method to utilize 
for exploring the motion characteristics of three degree 
of freedom airframe-autopilot systems, and was the first 
approach to airframe-autopilot system analysis.* ° Its 
utility stems from the fact that all automatic control 
systems involve control deflections which are made 
functions of airframe motions. 


ative’”’ approach. 


These control deflec 
tions impose forces and/or moments upon the craft 
which are proportional to airframe motion quantities. 
The action of the control system in producing such 
forces and moments can then be thought of as modi- 
When 


automatic control systems which effectively create new 


fying or creating aircraft stability derivatives. 


derivatives are considered, the equivalent-derivative 
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approach tends to be less useful. One can, of course, 
derive new approximate factors for different and higher 
order ‘equivalent airframes,’’ but the labor becomes 
excessive. Also for situations where the dynamics of 
the control-sensor combination greatly modify the 
simple proportional relationship between the sensed 
motion and the output control deflection, the equiv- 
alent-derivative approach becomes invalid. In such 
instances, servoanalysis techniques are eminently more 


suitable, and the application of these methods will be 
treated later in this paper. 

As an example of the equivalent-derivative approach, 
consider a control system consisting of a rate gyro, 
properly oriented to measure the yawing velocity, 
feeding a signal to servo actuators which move the 
rudder in response to the gyro signal. If the servo 
and sensor dynamics are neglected and the system 
assumed linear, the controller transfer function is 


TABLE 5 
Lateral Numerator Approximate Factors—Aileron 
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Lateral Numerator Approximate Factors—Rudder 
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5,/r = K,. Then all three lateral equations of motion will be modified by the control system—viz., 


(s— Y,)6 — 
—N,'B — N,'p 


The total effect of the control system is to augment 
L,’ and N,’ and to modify the side-acceleration term 
which originally came from the centripetal force m Ur. 
In most practical cases, however, the principal effect 
of the r to 6, feedback is to augment N,’. Referring to 
Table 4, which gives the lateral-denominator approxi- 
mate factors, it is seen that an increase in — NV,’ by an 
increment, —AJN,’, will (1) increase the Dutch roll 
damping by the amount of — AN,’; i.e., 


A2(tw), =—AN,’ 


(2) make the spiral mode more stable (assuming that 
L,’'/Ng’ is negative) by adding an increment 


A1/T. = —Tp,(g/Us)(Le'/Ne') AN,’ 


These changes to the denominator factors are equivalent 
to those which would be found by closing the control 
system loop using servoanalysis methods. Therefore, 
the closed r to 6, system is characterized by an increased 
Dutch roll damping and a more stable spiral mode when 
the sign of K, is properly taken. The converse will 
occur for an opposite sign of K,. 

Besides the changes to the poles of the lateral transfer 
functions given above, an increase in —JN,’ will also 
modify many transfer function zeros. For instance, 
from Tables 5 and 6 it can be seen that the numerators 
of the control deflection transfer functions are affected 
as follows: (1) the roll-numerator damping, 2¢,w,, will 
be increased for the aileron case by that increment A/V, 
of AN,’; i.e., 


AK yatdyg = —AN, 
and for the rudder case by the increment AN,’; i.e., 
A2f,.w,, = —AN,’ 


(2) the damping of the quadratic numerator for N,(s) 
is increased almost directly by an incremental change 
in —AN,. 

All yaw-rate numerators are unaffected by changes 
in — J,’ since no subscript r derivatives can be present 
in the numerator determinant J,(s). 

The above example demonstrates the application of 
the approximate-factor formulas to the assessment of 
relatively low-gain feedback in an r to 6, control system. 
Similar examples could have covered p, g, 7, u, v, or w 
signals fed to appropriate surfaces. While the example 
given was for a single-sensor control loop, the equiva- 
lent-stability-derivative method is not restricted to such 
single feedbacks. Indeed, in practice, the method is 
probably most valuable when several feedbacks are 
used, since these system types are simply handled by 
the equivalent-stability-derivative procedure but be- 
come a complex multiloop analysis problem when servo 
techniques are employed. Applied in preliminary 


(g Use + 
—L,'8 += L,’)p — 
+ [s— (N,' + N,’K,)\r = 08 


(1 — Y,,*K,)r = 0) 
(L,’ + L;,’K,)r = 0> (20) 


analysis of multiloop problems, the practical utility of 
this method is sufficiently great to justify the develop- 
ment of airframe approximate factors for use in it 
alone. 


Closed-Loop Servo Implications of the Airframe 


Precise analysis of airframe-controller systems re- 
quires the use of servo methods; the optimization proc- 
ess requires consideration of the airframe as an alter- 
able element. In practice, the determination of the 
closed-loop effects due to variation in the airframe ele- 
ment demands a separate analysis for each of several 
sets of relative pole, zero locations consistent with vary- 
ing combinations of airframe and control systems. The 
results of these analyses provide the foundation re- 
quired to develop an appreciation for the potentialities 
of tailoring the airframe-controller system. The rela- 
tive roles to be played by each of the system elements 
can then be determined on a rational basis by consider- 
ing competing possibilities. 

The series of servoanalyses referred to above is 
straightforward as far as procedure is concerned, but 
complicated in detail. Therefore, in considering the 
servoanalysis approach, a ‘“‘single-sensor control loop 
system”’ will be used here to provide a simple example. 
The types of systems classified under this title have the 
distinguishing characteristic that their functional 
operation can be described by single-loop block diagrams 
composed of analytically simple transfer function 
blocks. Such a system is shown in Fig. 1. The for- 
ward-loop transfer function, Y,(s), may include closed- 
loop transfer functions of simple inner loops such as a 
positional servo. Also, the feedback transfer function, 
Y;(s), may be a combined representation of two or more 
sensors and additional simple closed-loop transfer 
functions. In both of these cases the transfer functions 
resulting from the sensor combination or the inner loop 
closure are simple, and the overall outer control loop, 
or “sensor loop,” is a ‘‘single control loop.’’ In other 
words, while each of the single series blocks in the block 
diagram of Fig. 1 may involve more than one loop 
and/or more than one sensor, the system is still a 
“‘single-sensor control-loop system.” 

Further restrictions are desirable because of the 
complex form of airframe transfer functions and the 
nature of closed-loop studies. In particular, both 
actuator and sensor dynamics are assumed to be 
ideal. These assumptions are made for two reasons: 
first, to simplify the analytical treatment, and second, 
to focus attention upon the airframe rather than the 
control system. In a practical flight-control design 
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problem such simplification, of course, would not be 
employed. 

In the diagram of Fig. 1, the two disturbance blocks 
essentially convert the disturbance input, y, into 
equivalent surface deflections appropriately ‘‘scaled”’ 
for ultimate conversion of yw into components of the 
airframe motion quantities § and. The total directly 
controlled airframe-motion quantity, A, is given by 


YA; J I Nw \ 
. = ee tS 
1+ ¥e¥A, “* ¥. My"f 
-| von |e 4 i | | (21) 
1+ You LY, Vat} Nys 
where Yo, = Y,Y;A; 


The components entering into the calculation of Af are 
seen to be (1) the “closed-loop’”’ transfer function, 
Yor/(1 + Yoxz) (this term contains all of the system 
(as opposed to input) dynamic and static character- 
istics), (2) the command input, \,, modified to system- 
output units by dividing by the feedback transfer 
function Y,, and (3) the disturbance, yu, with units 
changed to output terms by a suitable factor. 

The transfer functions Y; and Y,¥Yy are ordinarily 
relatively simple. In the discussion to follow they are 
further idealized. Also, the disturbance numerators, 
Ny,, are given in terms that are simplified, but suffi- 
ciently precise for most purposes, in reference 3, and 
the control numerators WN); have been presented in 
Tables 1 through 6 in various forms. Therefore the 
contribution of the terms listed as items 2 and 3 above 
is easily computed with the formulas already available. 
The closed-loop transfer function, on the other hand, 
has not yet been considered. Therefore the following 
discussion will be devoted to investigations of the 
influence of possible forms of specific airframe transfer 
functions upon the ciosed-loop transfer functions. The 
closed-loop characteristics will be simply and concisely 
shown in terms of the airframe transfer function by 
using standard methods of servoanalysis. 

When all of the component elements listed above are 
known, one can readily determine the response of the 
directly controlled airframe-motion quantity, A, to the 
various inputs. Further, the influence of possible air- 
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Fic. 1. Single-sensor control system with command and 
disturbance inputs. 


frame changes upon the component elements can be 
assessed through the medium of the airframe transfer 
function approximate factors. In actual practice, these 
pieces of knowledge allow the evolution of a rapid, cut- 
and-try, airframe-optimization procedure for single- 
sensor control-loop systems (as far as the directly con- 
trolled airframe-motion quantity is concerned). 

As the single-sensor control-loop example, consider a 
lateral control system in which roll angle and roll rate 
are both fed to the aileron. In terms of the block dia- 
gram of Fig. 1 and the expressions given in Eq. (21), 
this is equivalent to defining 


Y,=K,(1+7,s) and Y,=1 


Actually the basic roll and roll-rate system can incor- 
porate the 1 + 7\,s term in either Y, or Y;and systems 
have been built both ways. However, the location of 
the lead within the block diagram will be of no further 
concern here since the major interest is centered upon 
the closed-loop transfer function Yo,/(1 + Yorz). 


(1D) Closed-Loop Transfer Function Considerations Using 
Simplified Airframe Representation 

To introduce some of the concepts, it is advantageous 
to consider first a simplified roll transfer function, de- 
rivable from the basic equations of motion by setting 
the (s — Y,)8 term equal to zero. Reference 3 shows 
that such a simplification yields the roll-subsidence and 
spiral modes given in Table 4. For such a simplifica- 


tion, and letting Y;,* = 0, 


Liq’ (1 ax (Niq’ Le’ L,'Na’)| 


f _ { as r re? , re ’ 
ta = 52 I [Ny! — (g/Ue) \Lp’/No’) — Ly’) 5 + (/Us) IN;'(Lp’/Np’) — L,’] oe 
Kes, 


(s w))? + (2f1s ‘w) + l 


K gia 
Oo 
(T.s + 1)(Zes + 1) 


UsLag' No’ (1 —(Nig'Lp'/Ltg'Np’)} 
gN,’ Ly’ [1 — (Ng’L,’/Le'N,)] 


where Kew = 


+ A corresponding development for the indirectly-controlled 
airframe-motion quantities, such as & in Fig. 1, leads to the same 


implications. 


(t T,) + (1 Tr) = 2610 = 
—L,’ - in, = (g/Uo)| (Le’ Nz’) 


1/T.TR = wr’ = (g/U0)(Nr'(Lp'/N a’) — Ly’) 
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Since the case with two real roots is the more common, 
it will be used here. 

The open-loop transfer function is formed by taking 
the product of the airplane and controller transfer 
functions—viz., 

ae KK ei4(T 9s + 1) (23) 

Oo” (Ts + 1)(Trs + 1) 2 

Since the roll subsidence is always greater than the 
spiral root, the simplified system characteristics will be 
strongly dependent upon the relative location of 1/7), 
and 1/7. Fig. 2 shows generic open-loop Y(—<a) 
Bode diagrams® for the two possible cases (excluding 
the trivial one when 1/7, = 1/T,). Root locus plots’ 
corresponding to the Bodes are shown in Fig. 3. Al- 
though Y(jw) Bode diagrams are also useful representa- 
tions in the analysis procedure, they add little to the 
present discussion. 

Each one of the three diagrams mentioned contains 
all the information available concerning the control 
system and its open- and closed-loop transfer functions 

i.e., the closed-loop transfer function poles and zeros 
can be obtained from, or are presented on, each of the 
plots. The reason for using al three forms in normal 
generic analyses is that each plot emphasizes certain 
features of interest more than the others. The plots 
are supplementary in that they each yield the closed- 
loop secrets with various degrees of difficulty in various 
regions. While it is not within the scope of this paper 
to describe the use and interpretation of these plots in 
detail, the two most applicable to the present situation 
will be discussed rather fully for illustrative purposes. 

(a) The root locus plot, Fig. 3, consists of loci which 
are paths defining all possible closed-loop transfer 
function poles. The plot uses a linear scale, and the 
roots are shown as complex numbers with the abscissa 
giving the real part and the ordinate the imaginary part 
of a closed-loop root. Gain is a parameter on the loci 
(the arrows are drawn in the direction of increasing 
open-loop gain). 

(1) At zero gain, the poles (X’s) of the closed-loop 
transfer function are simply the poles of the airframe 
transfer function alone. In other words, with the 
control system gain turned down to zero, the airframe 
motions are those of the airframe itself. In Fig. 3, 
these airframe dynamics are defined by the spiral 
root (—1/7,) and by the roll subsidence root (—1/7). 

(2) At very latge controller gain, the roots defining 
the airframe’s motions fall along the high-gain, straight- 
line asymptotes to the locus, or tend to approach the 
zeros (©) of the open-loop transfer function. For the 
simplified situation of Fig. 3(a) there is one transfer 
function zero, and the high-gain asymptote, @), is 
actually part of the locus throughout its entire extent. 
This portion of the locus, and the other real root pro- 
ceeding along the path @), result from a decoupling at 
the encircled 6 point of the coupled spiral and roll- 
subsidence modes. 

(3) At intermediate values of gain, the closed-loop 
roots defining the airframe motions are somewhere on 











= <1 


Fic. 3. Root locus for simplified roll case with lead equalization. 


the locus between the airframe-alone poles and the zero 
and asymptote. In Fig. 3(a), the “‘effective’’ closed- 
loop spiral root becomes a larger negative number 
(path @) and the modified roll-subsidence root becomes 
smaller (path @)) until the two couple (at the encircled 
a point) to form an oscillatory pair. The locus of roots 
between the encircled points a and 0 is a circle. If 
the imaginary axis is y and the real axis x, the equa- 
tions for the locus can be shown to be‘ 


[x + (1/T,)]? + y? = (1/T,)? + (1/T eT.) 
(1/T,)((1/Tr) + (1/T.)] (24) 


y=0 (25) 


The circle described by Eq. (24) has its center at 
(—1/T7,, 0) and has a radius which is approximately 
V1/T, V1/T, — 1/Tpe (if 1/Te > 1/T.). When 
1/T., < 1/Tp, the radius is imaginary, and the situation 
corresponds to Fig. 3(b). Eq. (25) has regions of valid- 
ity regardless of the relative orientation of 1/7), 1/Tp, 
and 1/7,; it is therefore a companion to Eq. (24) when 
the circle exists, or the sole locus equation if the radius 
is imaginary. 

(b) The amplitude ratio portion of the Y(—a) Bode 
plot of Fig. 2 is obtained by (1) substituting the real 
variable, —o, for the complex variable, s, in the open- 
loop transfer function; and then (2) plotting 20 log 
|Y(—o)| versus log oc. The amplitude-ratio asymp- 
totes are identical to those of a Y(jw) Bode diagram. 
The phase angle is either zero or — 180 degrees, depend- 
ing upon the sign of Y(—c). The portions of the plot 
where the phase angle is —180° satisfy the angle 
criterion for a closed-loop root—i.e., for a closed-loop 
root to exist, V(s) + 1 must equal zero, so Y(s) = —1 
or ZY(s) = —180°. The amplitude ratio in the 
regions where the phase is —180° is therefore a locus 
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Fic. 4. Root loci of 66 = Ky(Tys + 1)gsystem for Condition I. 


bil. < 1/Te < 1/Tr. 


of the negatives of real roots versus gain. The Y(—<) 
Bode is both an open-loop plot and a locus of all real 
closed-loop roots. 

(1) When 1/7, > 1/7 [Fig. (2a) ], as gain is increased 
the modified spiral and roll-subsidence proceed down 
paths @) and (@) until the zero dd line reaches the en- 
circled a point. At higher gains the roots become 
complex. At still higher gains the presence of 1/7), 
becomes a more emphatic factor in the character of the 
system. For the zero dd line at the encircled } point 
or lower—i.e., higher gains—the two oscillatory roots 
have returned to the real axis: one goes into (—1/7,) 
while the other increases (negatively) with increased 
gain. 

(2) In Fig. 2(b), the presence of the zero, —1/7,, 
between the —1/7, and —1/7> poles leads to a system 
which exhibits considerable differences from the above. 
The behavior of the closed-loop roots as gain is increased 
is straightforward—i.e., the modified spiral goes down 
path @ to the lead term 1/7 ‘» and the modified roll 
subsidence proceeds along path (@). 

Paths @ and @) and encircled points a and 5} are 
correlated directly with the same paths and points on 
the root loci of Fig. 3. As can be appreciated from the 
above discussion, the airframe approximate factors are 
as useful in the servo techniques as they are for the 





equivalent-stability-derivative approach. In the root 
locus, for instance, they define the starting and end 
points of the locus; and the intermediate points are 
completely determined by a weighted effect of each 
pole and zero. For example, in the case discussed here 
(with 1/7, > 1/T7,), it is apparent that the response 
time, total system damping, etc., are direct functions of 
1/Tr. Since 1/7, is known approximately in terms 
of airframe parameters, the influence of these upon the 
closed-loop system can be readily assessed. 


(2) Closed-Loop Transfer Function Considerations With 
the Complete Airframe Representation 


If the complete three degree of freedom airframe 
characteristics are used, the ®;, transfer function, with 
Y;,,* = 0, has the form 

— K gal (5, Wy,)? + (25 oo /Wq)S - 1] 

(Tis + 1)(Trs + 1)[(s/wa)? + (2648/wa) + 1] 
(26) 


®;, = 


The approximate values of 73, 7p, and Ag,, for this 
complete case are the same as those previously used in 
the simplified roll-control situation. The primary 
change in going from the simplified case to the more 
complex one is, then, a multiplication of the simplified 
airframe transfer function by a ratio of quadratics. 

As noted in reference 3, the values of ¢(gwg and 
fe.%e, are ordinarily about the same; the ratio 
of w,, and w, can be either greater or less than unity; 
and 1/7, may be either greater or less than wg (or 
w,,)- Most of these relative pole and zcro locations 
are amenable to fairly simple changes by modifying the 
airframe parameters. The various possible combina- 
tions and permutations of these parameters, as well as 
the value of 1/7.,, are therefore of interest. 

As one might suspect from the discussion of Figs. 2 
and 3, the relative location of 1/7, 1/7,, and 1/7, 
can be used as a principal classifying device. Roman 
numerals will be used for this basic division—.e., 


Ve ee ty mee te 
1/T, < 1/Tp < 1/T, 


Condition I 
Condition II 


The practical single-sensor control-loop cases for 
Condition I, for all types of w,,/wg ratios, are shown in 
the root loci of Fig. 4. As gain is increased, the spiral 
root always goes toward the lead-equalization zero, 
the roll-subsidence root goes out the-real axis, and the 
Dutch roll roots move into the roll-numerator quadratic 
zeros. In the latter portion of the locus, the modified 
Dutch roll roots always approach the roll-numerator 
quadratic along counterclockwise paths which resemble 
a sector of a circle. It is clear from this behavior that 
situations corresponding to w,,/wg > 1 can become 
unstable (as in Fig. 4, curve a), while the converse is 
true when w,,/wg < 1. For the modified Dutch roll 
to actually become unstable, a phase dip in a Y(jw) 
Bode diagram exceeding —180° is required. The 
presence of such a dip depends, of course, upon the 
actual relative values of £4, £44) gq @a, 1/Tp, and 1/7). 
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Condition II situations can be much more varied 
than those of Condition I. The chief reason for this 
is the mutual influence and interaction of the coupled 
roll-subsidence and spiral mode with the w,,, ¢,, and 
wa, 4 pair. In the discussion of the simplified case, it 
was noted that the coupled spiral and roll-subsidence 
mode was a circle. The addition of the ratio of quad- 
ratics to the airframe transfer function modifies this 
somewhat, but a large part of the root locus is still 
quasi-circular. Now if the radius, y, of the circle gener- 
ated by the real poles and zeros (—1/7,, —1/Tk, 
—1/T,) is about the same as the distance, v, from 
—1/T, to the wg zero (if w,,/w¢ > 1) or the w,, pole 
(if w,, ‘wa < 1), then one would expect strong inter- 
actions between the quasi-circle due to —1/7,, etc., 
and the loci which in Condition I departed from the 
Dutch roll roots and arrived at the numerator quad- 
ratic zeros. Just such interactions do exist; and so 
the relative magnitudes of 7 and v become an appro- 
priate classifying characteristic. 

The distance v, from —1/7., to a point —fw,, 
wo,W1 — ¢*? (where w, and ¢ can refer to either the 
@,, Zeros Or wg poles, whichever is physically closer to 
1/T,,) is 


v2 = [(1/T,) — Sen]? + (wn V1 — ¢2)) 





. 7 
(1/T,)? — 2etn/Tp + wn? ; oo 

The approximate radius of the circle, 7, is 
n? = (1/T,)? — (1/T,T pe) (28) 


By combining Eqs. (27) and (28), the classifying cri- 
teria become 
n° 5 v?, or (1, Tn) [C1 ‘T R®n) ~ 2 ] = 
m<v’, or (1/Ton)[(1/Tron) — 26] > 
Besides the “‘influence’’ parameter discussed above, 


the w,,/wg ratio is also of consequence. Therefore, the 
total cases of some practical interest for Condition IT 


include 
wy,/wa <1 wy,/wWa > 1 


Case [:.. n< » Case 1’: n< pv 
Case 2: n= py Case 2’: 9 = p 





Case 3: >vp Case 3’: n>v 

The unprimed cases (w,,/wg < 1) are shown in Fig. 5. 
In Case 1, the circle centered approximately on —1/7), 
has little influence on the shape of the other loci. In 
essence this case is similar to that shown in Fig. 3 with 
the addition of a quadratic pair which are widely 
separated from the other roots. 

Case 2 is a transition configuration. The circle is 
still largely present, but is broken up in the neighbor- 
hood of the w,,, w¢ pair. The coupled roll-subsidence 
and spiral mode starts along the circle, and then takes 
off in a different direction toward the roll-numerator 
zeros. The modified Dutch roll roots, as one proceeds 
along the locus from low to higher gains, seem to be 
starting toward w,, as in Case 1, and then change direc- 
tion to go into orbit about —1/7,,. 

Case 3 shows the situation after transition. The 
Dutch roll is now the basis of the circular segment, and 
the coupled spiral and roll-subsidence mode advances 
directly into the roll numerator. 

The primed situations ((wy,/w¢ > 1) are surveyed in 
Fig. 6. An example of Case 1’ is not shown since it 
would be similar to Case 1 except that the locus starting 
at wq and moving toward wy,, while still counterclock- 
wise, would be rotated roughly 180° (since wy, > wa) 
thus bringing portions of it into the right half (unstable) 
plane. 

The example shown for Case 2’ is not as interesting 
as its counterpart, Case 2, since no apparent transition 
effects are evident. The locus shown does provide a 
good illustration of how close one can come to a transi- 
tional situation with little apparent effect upon the 
locus. 

The Case 3’ example is similar to that for the un- 
primed case with the notable addition that instability 
of the coupled roll-subsidence and spiral mode can 
conceivably exist. In these primed cases, as with the 
previous instances where w,,/w¢ > 1, the simplest 
quantitative picture of the stability problem is given by 
the phase dip on the Y(jw) Bode representation. 

For all the Condition I and II cases discussed, the 
high-gain behavior of the closed-loop transfer function 
will be essentially the same—i.e., 


T . 2 9¢ s/w 
(7's + 1) [(s/wy,) + (25 645 W oq) + 1] (30) 





= You ——" | =e | ? 
lL + Fe 1+ KyKos.) (Tors + 1)(Tgep8 + 1) [(s/we,)? + (2Fers/wer) + 1] 


which, as A,K,;, becomes much larger than unity, 
approaches 


[Yor/(1 + Yor) Ik Ke, targe = 1/(Texrs + 1) (31) 


since, then, 7,,, = Ty, wer = wg. and fcr = fy, 


The value of 1/7¢, at high gain is simply the inter- 
section of the zero db line and the high-frequency asymp- 
tote of the Y(jw) or Y(—o) Bode diagrams. 

The cases presented above as examples are also useful 
for other purposes. For instance, an important ob- 


jective in the design of a simple and reliable control 
system is to avoid changes in the equalization, 1/7), 
throughout the flight regime. As the flight condition 
is changed, the airplane characteristics also change; 
as a result, the relative location of the controller-air- 
frame system poles and zeros may be altered drasti- 
cally. In Fig. 5, for instance, Case 1 could be a high- 
speed condition, Case 2 could represent intermediate 
or cruise conditions, and Case 3 a post-launch or low- 
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Fic. 6. Root loci of 66 = Ky(Tys + 1)¢ system for Condition II (1/7, < 1/Tx < 1/T,) when we,/wa > | 
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speed condition, with the control system 1/7, held 
constant. As can be seen from this sequence, the 
character of the closed-loop response could change 
drastically as a function of the actual flight condition. 


(3) Summary of Significant Closed-Loop Considerations, 
e(T,s + 1) > & 


Pertinent points developed in the text above are 
summarized below. One should realize that 1/7,, 
fa a $4, and w,, are airframe parameters that are 
relatively variable (although not necessarily inde- 
pendently variable). 

(a) Spiral, 1/T7,—The spiral mode has little effect 
upon the lead-equalized roll system studied for phys- 
ically reasonable gain values, although a divergent 
spiral places a lower limit on suitable gain values if it 
is to be stabilized. 

(b) The Ratio of Lead Equalization to Roll Subsidence, 
T,/Tz—The ratio of the lead-equalization time con- 
stant, 7, to that for the roll-subsidence, 7’, is a suit- 
able classification parameter for lead-equalized roll 
control systems. 

For 7/7 < 1, (a) all systems examined have real- 
axis root loci in which —1/7, goes to —1/7\,, and 
—1/T,p advances toward —@© as K, increases (the 
roll-subsidence and spiral modes do not couple), and 
(b) the Dutch roll roots proceed in a counterclockwise 
direction into the roll-numerator zeros. 

For 7,/Tr > 1, (a) the roll-subsidence and spiral 
modes couple and break off the real axis, sometimes in a 
quasi-circular path back to the real axis and thence to 
—o and —1/7.,, and sometimes to the roll-numerator 
zeros, and (b) the Dutch roll roots either go into the 
roll-numerator zeros along a counterclockwise path or 
form a quasi-circular segment to the real axis (and 
then to— and —1/7,). 

(c) The (1/Ton)((1/TRwn) — 2¢] Influence Param- 
eter—This ‘“‘influence’’ parameter is of some conse- 
quence when 7\,/7p > 1. The w, and ¢ refer to either 
w,, and ¢,, or wa and ¢4, whichever set represents the 
point closest to 1/7.,. The relative magnitude of the 
parameter compared with unity offers a crude and 
approximate classification device to separate the de- 
tailed locus interaction forms appearing at intermediate 
gains. When 7/7, > 1, and 
(1) (1/T on) [(1/T rn) — 2] <1, (i.e., 9? > v?) 
the Dutch roll roots will ordinarily go toward the real 
axis along a quasi-circular segment; the coupled spiral 
and roll-subsidence roots proceed to the roll-numerator 
zeros; 

(2) (1/T yon) [C/T ron) — 26] > 1, (ie., 9? < v?) 


the Dutch roll roots will ordinarily move to the roll- 
numerator zeros and the coupled spiral and roll-sub- 
sidence will form a quasi-circular locus. 


(d) Zhe Ratio of Roll-Numerator Frequency to Dutch 
Roll Frequency, w 


‘wa—A ratio of w,,/wa¢ less than 


Ya 


unity is preferable in all cases and mandatory in some 
to retain stability at intermediate gains. 


Examples of Competing Systems 


The various forms of the controller-airframe transfer 
function given in the foregoing discussion indicate the 
possibility of setting up a number of competing systems. 
This possibility can best be demonstrated with a hy- 
pothetical example considering an airframe which 
initially has characteristics with certain common defi- 
ciencies from a single-sensor control viewpoint. By 
relating these shortcomings to the airframe configura- 
tion, certain changes in geometry can be investigated 
to determine their influence on the required controller. 
The competing systems can then be compared to deter- 
mine possible advantages. 

Consider the roll control of an airframe having very 
lightly damped Dutch roll and numerator modes, and 
having all other dynamic characteristics quite similar 
to those shown for Condition Ia of Fig. 4. Such a 
situation frequently arises in practical cases for flight 
conditions corresponding to high-altitude cruise oper- 


‘ations. Since the assumed poor Dutch roll damping is 


lower than that shown in Fig. 4 (which corresponds to a 
¢ = 0.1), the root locus connecting the Dutch roll pole 
with the numerator zero will pass into the right half 
plane (negative damping) at fairly low open-loop gains. 
Closure of the roll loop at any gain whatever will de- 
crease the damping of the Dutch roll mode which is 
already too low. 

A conventional solution to such a problem is to add 
a yaw-rate-to-rudder loop to the controller, which, so 
far, consists of only the roll loop. With this auxiliary 
control, certain of the ‘‘augmented airframe’’ poles and 
zeros will shift to the left, as already noted in the dis- 
cussion of the equivalent-derivative approach. The 
important displacements in this case are those for fqwa 
and {,w, (respectively, the Dutch roll and numerator 
dampings) which, for rudders with low rolling moments, 
will be about equal. If these roughly equal shifts are 
made sufficiently large, the resulting root locus loop 
shown for Condition Ia of Fig. 4 will always be in the 
left half plane. The complete closed-loop Dutch roll 
will then always be positively damped, even though 
the presence of the roll-numerator zero will reduce this 
damping somewhat depending on the roll gain em- 
ployed. 

As competition for the foregoing solution, it is per- 
tinent to investigate the airframe configuration changes 
required to achieve a pole, zero relationship corre- 
sponding roughly to Condition Id of Fig. 4. If this 
could be simply accomplished, then as indicated by the 
root locus plot, the system would be stable for all roll- 
loop gains without requiring the auxiliary rudder loop. 
Furthermore, the Dutch roll damping could be con- 
siderably increased by proper choice of gain. The 
basic change required to achieve the desired approxi- 
mate pole, zero configuration involves reducing, to less 
than unity, the ratio of 

















W 
inc 
fer 
ne: 
Sir 


1 
po: 
ori 
ail 
ail 
inc 


giv 
ac 


hec 
sin 
eff 
ma 
the 
the 
car 
ine 


Cn; 


vali 
crea 
can 


(Ta 


erro 





ome 


sfer 
the 
ms. 
hy- 
ich 
efi- 
By 
ira- 
ted 
ler. 


‘s- 














SO 


— 


OPTIMIZATION OF THE 


(27) #1 - Te 
Wad | F Se ’Ns 4 
oe [(Cns,/ Cts.) + (r2/Iz)|[(Cig/Cng) + (22/12) ] 
[1 + (2 T,)(Cig Cug) | 
(32) 


In the root locus plot of Fig. 4, the value of (w,/wa)? 
corresponding to Condition (a), the original configu- 
ration, is about 1.25, whereas a more desirable value, 
corresponding to Condition (d), is approximately 0.6. 
In setting up the airframe characteristics for the original 
example, factors which contributed strongly to the re- 
sulting poorly damped Dutch roll were a highly effective 
dihedral (large negative Cig) and a moderately low 
(but negative) value of J,,. The influence of these 
parameters on the Dutch roll damping can be appre- 
ciated by converting the approximate expression of 
2¢awa given in Table 4 to basic terms—viz., 


—_ pSl ly J z ] (*)( oe ) — 
2m (os + 2\k, Cn + 2 o 


(i. T,)(Cig Cng) + Ce Fr.) 
1 + (J,2/I,)(Cig/Cng) 


l b P | 
E > x) (c., +; c,)|t (33) 


With Cug positive (which is normally the case and is 
indicated by the moderately large value of wg) and re- 
ferring to Eq. (32), the last numerator bracket will be 
bracket positive. 


2Fawa = 


negative and the denominator 
Since w,/wa > 1 for this example, it follows that the 
term (Cn;,/Cis, + I,:/],) must be positive. Since 
I,./1, is negative, this means that Cn;,/Ci;, must be 
positive and larger numerically than /,./J,. The 
original configuration is thus seen to involve favorable 
aileron yawing moments (e.g., corresponding to inboard 
ailerons),} high effective dihedral, and a low positive 
inclination of the principal axis (low negative /,.) for the 
given flight situation. 

To reduce w,/wq to the desired value it is clear that 
a change in sign must be effected for the term, Cn;,/ Cis, 
+ I,./I, Decreasing, but maintaining positive, di- 
hedral would not accomplish the desired reduction 
since, as Cig — 0, w,/wa —~ 1. Also a reduction in 
effective dihedral may not be actually possible since a 
major portion of Cig may be due to effects other than 
those associated with geometric dihedral. There are 
then two possible basic configuration changes which 
can be made. These are associated with either (1) in- 
increasing the negative value of /,,, or (2) making 


Cn;, more negative. The first involves raising the 


t Inboard ailerons tend to violate one of the conditions of 
validity required for Eq. (32) since such an arrangement usually 
creates Y5q* terms due to sidewash effects. The errors involved 
can be checked in detail on the basis of the rudder approximation 
(Table 6) which includes all effectiveness terms (Ls, Ns, and 
Y5*), In the specific numerical case later considered these 
errors are small, so Eq. (32) is retained for simplicity. 


FLIGHT-CONTROL, AIRFRAME SYSTEM 


to 
NJ 


TABLE 7 
Variation of (wy/w,)* With Stability Parameters for a 
Hypothetical Case 


T ~~ 4 \ 


. g xz ‘ “sg 
iguratior —=9; rt = - 53 tT “(1 - hh + 
\ x 


principal axis of inertia; the second involves either 
moving ailerons outboard or coupling them to the 
rudder. 

To give a better appreciation for the magnitude of 
the changes involved, consider the numerical values 
given in Table 7. As indicated by the table, the de- 
sired change in the character of the root locus for the 
hypothetical example can be accomplished by either 
(1) increasing the principal axis flight attitude by 0.05 
radian (about 3°) (a) through a decrease in wing 
incidence, or (b) through a mass redistribution effort; 
or (2) decreasing the favorable aileron yaw to approxi- 
mately zero (or less) by (a) shifting ailerons outboard, 
or (b) cross-feeding aileron signal to rudder. 

The elements for optimizing the airframe and roll- 
control system are now present in that a total of five 
competing arrangements have been set up as possible 
solutions to the basic problem. These all involve the 
basic equalized-roll-to-aileron loop and (1) the original 
airframe plus a yaw-rate-to-rudder loop, (2) the 
original airframe plus an aileron cross-feed to rudder, 
(3) the original airframe with modified mass distribu- 
tion to achieve A/,,/J, = —0.40, (4) the original air- 
frame with wing incidence reduced approximately 3° 
to achieve AJ,,/J, = —0.40, and (5) the original air- 
frame with ailerons shifted to an outboard position. 

Since, in this crude example, no attention has been 
paid to considerations other than the improvement of 
closed-loop Dutch roll damping, the competing sys- 
tems may not be equivalent when all modes of dynamic 
performance are considered. This reservation applies 
particularly to the possible adverse effects of sensor and 
actuator dynamics, which have been completely neg- 
lected, and to various steady-state problems. Despite 
these reservations, the competition appears to be valid 
from a system-synthesis standpoint. The basic choice 
of one or the other of the foregoing solutions must then 
depend on the optimization criteria employed and the 
physical implications of the various possibilities in 
regard to performance, weight, reliability, cost, etc. 
For example, if weight, reliability, complexity, main- 
tainability, or cost were the most important optimiz- 
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ation criteria for the example noted, the two competing 
systems requiring a rudder axis of control would almost 
certainly be put into an unfavorable position relative 
to the three possibilities where control by aileron alone 
could be sufficient. 


Preliminary Design Implications 


The case outlined above is, of course, intended to 
serve only as an example of the advantages acquired in 
problems concerning the airframe-controller system 
simply by considering possible airframe alterations. 
Many other considerations leading to potential system 
simplification and optimization can be deduced by 
-arrying out procedures similar to that used above 
for other system types, and then combining the an- 
alytical operations with detailed examination of the 
airframe transfer functions. For example, possibilities 
for tailoring the airframe element include: 

(1) Airframe basic geometry changes to modify such 
parameters as N,’, Ly’, and M,, thereby influencing 
controller static and dynamic requirements as well as 
affecting numerator characteristics such as the w,/w, 
ratio discussed previously. 

(2) Selection of appropriate  control-application 
points and magnitudes—.e., the various control arms 
(both positive and negative) and effectiveness—to 
obtain favorable airframe transfer characteristics from 
the standpoint of control. 

(3) Attempts to attain favorable airframe mass- 
distribution characteristics, particularly in /,,, which 
has a profound effect on both numerator and de- 
nominator terms. 


MARCH, 1960 


(4) Selection of sensor orientation and _ location 
within the structure to obtain favorable effects that 
would otherwise require complex controller-equalization 
procedures—e.g., use of tilt in lateral gyros—or place- 
ment of accelerometers near the center of percussion. 

Full utilization of these and other possibilities requires 
the basic understanding sought, and hopefully ad- 
vanced, by the work reported here. In addition it re- 
quires concerted, systematic teamwork between the 
stability and control aerodynamicist and the flight 
control designer at the earliest possible stage of pre- 
liminary design. 
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Heat Conduction in a Melting Slab' 


STEPHEN J. CITRON* 


Columbia Unwersity 


Summary 


The behavior, after the initiation of melting, of a slab insulated 
on one side and subjected to a heat input varying arbitrarily 
with time on the other is studied. A method of successive ap- 
proximations is developed which requires the solution of an 
ordinary nonlinear differential equation; this equation is readily 
solved numerically on a desk calculator. Since any heat input 
is permitted, solutions by this method may be used to match 
conditions resulting from the aerodynamic heating problem. 
Curves of melt-thickness versus time are presented for a constant 
heat input and jor a typical heat input incurred by a body 
descending through the atmosphere. 


Symbols 


\ = space variable 


t = time variable 

k = thermal conductivity 

( = specific heat 

p = density 

K = k/pc, thermal diffusivity 

3 = latent heat of melting 

l = slab length 

Q(t) = heat flux 

Oo = reference value of the heat flux 

7,(x, t) = premelting slab temperature distribution 

T(z, t) = slab temperature distribution after melting has 
begun 

tia = melting temperature 

e = time at which melting begins 

ty = time for slab to melt completely 

s(t) = position of solid boundary after melting has begun 
(amount of material melted ) 

S = s/l, nondimensional form ot s 

3 = (x — s)//, nondimensional space variable 

T = x(t — t*)/l?, nondimensional time variable 

Tl = nondimensional time after melting has begun 
at which slab is completely melted 

(Tie = approximate value of 7; obtained from solution 

6 = 7 /T*, nondimensional temperature variable 

2M/r'/? = cT*/L, ratio of heat content at melting temper- 
ature to latent heat of melting 

r = kT*/iQ, ratio of heat flux k7T*/] to reference 
value of the heat input flux Qo. 

4 = 0/0: 
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Introduction 


rT RECENT YEARS, since the advent of hypersonic 
missile flight, the so-called ‘“‘re-entry’’ problem has 
been given an increasing amount of study, among the 
most important aspects of which has been that con- 
cerned with the reduction of structural damage caused 
by aerodynamic heating. To accomplish this pro- 
tection a variety of design methods have been discussed 
in the literature and used in practice; these include, for 
example, the use of heat sinks,' of transpiration cool- 
ing” and of mass-transfer cooling.* An accurate anal- 
ysis of the latter two methods would treat the inter- 
action of a gaseous region, a liquid region, and a finite 
solid region. Since the rate at which the liquid is 
formed is not necessarily equal to the rate at which it 
is swept away by the flowing gas, the dimensions of 
these regions will vary with time and, thus, ablation 
of the melt formed is an important part of the total 
problem. The question of ablation will occupy our 
attention in this investigation. 

The classical work in heat conduction contains many 
references! to moving (free) boundary problems, 
often called Stefan problems. It should be noticed 
that in all these examples the overall dimensions of the 
material under consideration remained constant while 
a phase-change front moved through the material— 
that is, no ablation was assumed to occur. In addition, 
in the solutions of the above problems the methods 
used do not permit arbitrary initial temperature dis- 
tributions prior to melting; this must, however, be 
allowed in the ablation problem to account for the heat- 
ing history before melting begins. 

The earliest work relating to ablation which over- 
comes these difficulties is that of Landau’! who studied 
in detail the case of a melting or sublimating semi- 
infinite solid under constant heat input. Under the 
assumption of total instantaneous removal of the liquid 
or gas formed, curves of thickness melted and melting 
rate versus time were obtained numerically by digital 
computer over a complete range of physical parameters. 
The assumption of total ablation implies that the pro- 
tective effect of the actual liquid region adjoining the 
solid body will be small—i.e., the heat input to the 
body will not have been reduced in passing through the 
liquid. 

The conditions under which this assumption is valid 
have been studied by Goodman’? who discusses ana- 
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Fic. 1. Slab geometry. 


lytically, under certain restrictions, the ablation due to 
aerodynamic heating. The major restrictions intro- 
duced are that the liquid have a Prandtl number which 
is small compared to unity, that the ratio of melt thick- 
ness to boundary-layer thickness be small compared to 
unity, that the solid region be semi-infinite in extent 
and that its temperature be held constant at its melting 
value. 

In reference 13, Sutton solves a similar problem, 
by use of a REAC differential analyzer, without some 
of the restrictions which were required in the analytical 
work of Goodman; the solid is, however, still semi- 
infinite in extent. A solution is obtained for a glass of 
high Prandtl number and the results are compared 
with those of Landau. For the particular example 
solved it is shown that only 70 per cent of the energy 
arriving at the gas-liquid interface reaches the solid 
body. Conversely, it is indicated that, for liquids 
with low Prandtl numbers, the shielding effect of the 
molten film is much smaller and may therefore be 
usually neglected in general. In the latter case the 
assumption made by Landau of instantaneous ablation 
is thus justified. 

As has been mentioned, in the work on ablation of 
Landau and the related work of Goodman and Sutton 
it was assumed that the solid material was semi-infinite 
in extent. This assumption, however, will lead to 
melting rates that are too slow—i.e., not conservative— 
because of the finite dimensions of the actual structure. 
The present investigation concerns itself with heat con- 
duction in a melting slab. The behavior, after the 
initiation of melting, of a slab insulated on one side and 
subjected to a heat input varying arbitrarily with time 
on the other is studied. A method of successive ap- 
proximations is developed and solutions are obtained 
through use of a desk calculator which permit an arbi- 
trary initial temperature distribution (to account for 
the heating history prior to melting) and an arbitrary 
heat input (to match conditions arising from the aero- 
dynamic problem). Solutions for a constant heat in- 
put and for a typical heat input incurred by a body 
descending through the atmosphere are obtained. 


Formulation and General Properties of the 
Solution 


Consider the slab of thickness / shown in Fig. 1, in- 
sulated at x = / and subjected to a net heat input Q(‘) 
atx = 0. We will let ¢ = ¢* be the time at which the 
temperature at x = O reaches the melting value 7%, 
while the interior temperature at this time will be 
To(x, t*). T(x, t) denotes the solution to the heat con- 
duction problem in the slab prior to melting (0 < ¢ < /*) 
obtained from the heat conduction equation: 


k O°T)(x, t)/Ox® = OTy(x, t)/Ot, O<x<l (1) 


under the conditions 


( alr kOT\| 

a) = = 

, Ox x=0 
Ox rt | 

(c) Ty(x, 0) = 0 


The first condition gives the heat input to the slab at 
x = 0, the second requires the slab to be insulated at 
x = 1, while the third condition prescribes the initial 
temperature distribution, taken here to be zero. 

Since one face of the slab is insulated the addition of 
more and more heat will eventually cause the slab to 
melt. After melting has begun, the position of the hot 
face is no longer fixed but will be an unknown function 
of time which must be determined and, hence, gives rise 
to a floating boundary value problem. The position 
of the solid boundary after melting has begun will be 
denoted by x = s(t) while the temperature at this 
boundary is held fixed at the melting temperature 7*. 

Mathematically the problem may be formulated in 
the following manner for the post-melting period, ¢> /*, 
se: 


x O?T (x, t)/Ox? = OT)(x, t)/dt, s(t)h<x<l (3) 


with conditions 


(a) Q(t) = —(ROT/Ox) + pLs\ 
(6) T(x, t) = T* f 
(c) Of (x, D Ox| 21 = 0 (4) 
(da) Te. t*) = Tole; 1°) 

(e) s(t®) = 0 


x = s(t) 


The first three are boundary conditions which pro- 
vide (a) that a part of the incoming heat at the melting 
boundary enter the solid while the remainder go toward 
overcoming the latent heat of melting L, (b) that the 
melting temperature exist at the melting face, and (c) 
that the slab be insulated at the other face. The last 
two are initial conditions which provide (d) that the 
temperature be continuous when melting begins, and 
(e) that the solid be of length / at this time. 

The presence of the moving boundary s(¢) thus intro- 
duces an additional unknown which necessitates three 
boundary conditions (4a, b, c) and two initial condi- 
tions (4d, e) in place of the two boundary conditions 
(2a, b) and one initial condition (2c) of the premelting 
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HEAT CONDUCTION 


problem. Certain general results can be immediately 
obtained from Eqs. (3) and (4a-e), and will now be 
discussed. 

(a) We see from (4a) that if Q(t) is to be continuous 
at ¢ = ¢* the melting rate s(t*) must, in general, vanish 

i.e., s(t*) = 0; L # 0. In the special case, however, 
where the latent heat is identically zero, a finite initial 
melting rate will not disturb the continuity of Q(¢) at 
t= t*. 

(b) As the surface at x = / is insulated the addition 
of heat to the slab will eventually cause it to melt com- 
pletely; the time (/,) at which this occurs may be de- 
rived from the set of Eqs. (3) and (4a-e) as follows. 
We proceed by integrating both sides of (4a) between 
‘* and ¢, with the result 


t t 
QO(8)deB = -{ k OT(s, B)/OxdB + pLs(t) 
ed {* t* 
The first term on the right may be rewritten using 
the fact that 
—k O1(s, t)/Ox = 
1 
(0 an] f pc] (x, t)dx + pc T* x | 


s(t) 
obtained by integrating Eq. (3) between s(t) and /, 
taking the derivative with respect to time outside the 
integral and using conditions (4b, c, e). 
One then has 


[ Q(8) dg = 


Jil 


f (0/08) | f pc T(x, B)dx + pc r*s(8) |B +pL s(t) 
* s(8) 
which upon integration yields 


I. Q(B) dB = p s(t)[L + cT*| + 


I 1 
i) pc T(x, t) dx — f pc T(x, t*) dx (5) 
s(t) 0 


As s — 1; i.e., complete melting of the slab is ap- 
proached; we obtain 
I 


*7 
| Q(8) d8 = ol [L + cT*] — { pc 7,(¢, t*) de (6) 
e” ( 


) 


from which, for any Q(t) the melting time may be com- 
puted. For the particular case in which Q(f) = Qp, a 
constant, Eq. (6) reduces to 


i, = ag + pl| [L + er") Qo} _ 
i 
(1 on) f pc To(f, t*) dé (7) 
0 
An expression for ‘* may be derived in a similar 


manner from a heat balance for the premelting period 
O0<t< ¢* resulting in 


i} Q(8) dB = f oc Ty(¢, t*) dé 
0 0 


7 This was first obtained by Landau!! using Gauss’ theorem 
on the heat conduction equation in the x, ¢ plane. 
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for arbitrary heat inputs and 


t* = (1 on) f pc 74(¢, t*) dé 
0 


for a constant heat input, Q. Hence, the total time 
required for a slab to melt completely under a constant 
heat input is given by 


t, = pli (L + cT*]/Qo} (8) 


(c) To ensure that the temperature in the interior 
nowhere exceeds 7* the temperature gradient at x = 
s(t) must be negative or at most zero, O7(s, t)/Ox 
0. From Eq. (4a) we see that this implies for a positive 


heat input the relation 
§ © Q(t)/pL (9) 


Since the gradient at the front face approaches zero as 
s — / we see that the equality holds at¢ = ¢,. Further, 
if Q(t) > O, the melting rate will have its maximum 
value at? = ¢,. 

(d) The solution of the postmelting problem for the 
limiting case J = 0; i.e., an infinite latent heat of 
melting; has not been presented here because of space 
limitations but could be obtained in a manner analogous 
to that of reference 11. For the other limiting case 
M = o~;i.e., no latent heat of melting; the initial melt- 
ing rate has been found to be not zero but 

S(0) =r 0*6 
027|,-0 
where the symbols are defined below. 

We will now present the problem in a nondimensional 
form. For the postmelting period it is convenient to 
change the independent variables from x and ¢ to z and 


T. 
z= [x — s(é)]/l, +r = («/P)(t — t*) (10a) 


We thus see that when x = s(t), = 0, and, hence, this 
transformation fixes the melting boundary. The 
transformation in time is used to shift the time scale 
so that it is zero when melting starts, therefore 7 > 0 
after melting has begun. The nondimensional temper- 
ature and the nondimensional melt parameter are de- 
fined to be 


6(z, r) = T(z, r)/T*, S(r) = s(r)/l (0b) 


respectively. The only physical parameters which will 
appear in the solution of the problem are the quantities 
M and r given here. 


M = (w'2/2)(cT*/L), r = kT*/IQ) (11) 


The parameter M may be interpreted physically as the 
ratio of the heat per unit volume required to raise the 
solid to the melting temperature to the latent heat of 
melting. The numerical constant 7'/*/2 in M was 
included for consistency with reference 11. The quan- 
tity Qo in r is a reference value of the heat input.f In 


t If Q(r)|..0 = O some other reference value may be taken 


without affecting the discussion. 
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what follows it will be convenient to take 


Qo = ae p=. ~ U 


Our problem then requires the determination of the 
temperature 6(z, 7) in the solid and the position of the 


moving free boundary S(r). 

In terms of these quantities, Eqs. (3) and (4a-e) be- 
come by direct substitution Eqs. (12) and (l3a-e), 
respectively, with corresponding physical interpreta- 
tions. 


2? + S$ d0(z, r)/dz = 
06(z, r)/O7r, O<z<1—S (12) 


070(z, 7), 


with conditions 


j 1/2 » 
(a) Q(r) cecunk 06 4 T ? $ 
Oo O2'. =0 2M 
(b) a0, 7) = 1 (13) 
(c) 06(1 — S, 7)/0z = 0 
(d) 6s, 0) = Tot, *)/T* 
(e) S(0) = 0 


The nondimensional equivalent of Eq. (9) is then 
S < (2M/w"/r)[Q(7r)/Qo| (14) 


We shall denote by 7, = (x//*)(t, — ¢t*) the dimen- 
sionless time after melting has begun at which the slab 
has been completely melted away. In the case where 
O(r) = Qo, a constant, for r> 0 we obtain 


l 
mr. = 7 t+ (9/?/2)(r7/M) -f Tle, 2")/T*)\ de (15) 
0 


by putting Eq. (7) in nondimensional form. 


Method of Solution 


To obtain the general solution it is useful to represent 
our nondimersional temperature @(z, 7) in the following 
form: 


6(z, r) = [To(z, ¢*)/T*] — 

[7y’(1 — S, t*)/T* ]e(z, 7) (16) 
where prime indicates 0/0z and ¢(z, 7) is now a new 
dependent variable. The reason for this change will 
become apparent as we proceed. Substituting Eq. (16) 
into Eqs. (12) and (13a-e) we obtain the equation and 
conditions which must be satisfied to determine ¢(z, 7) 
hence the temperature, and S(r), hence the amount of 
material melted. They are given below, each one cor- 
responding in physical interpretation to that of Eq. (3) 
and conditions (4a—e) which are listed in the same order 
and from which they stem: 





o*y § op (O/ Or)Ty’ | nek S, t*) a 

Oz? oz T)’(1 — S, t*) 
ov To" (2, t*) S73’ '(z, t*) (17) 
ye’ Ra-Eh * RA-sSPY 


and 


(a) Q(r)/Qo = (#'/?/2)(r/M)S — 

(r/T*)[To’(0, t*) — T’(1 — S, t*)¢’(0, 7)] 
(b) ¢(0, t) = O 
(c) dol — S, r)/Oz = 
(d) lim 7)’(1 — S, t*)¢(z, 7) = 0 


70 


(e) S(O) = O 


1 (18) 


Eqs. (17) and (18a-e) represent a nonlinear problem 
which we are unable to treat in an exact analytical 
manner. The alternatives exist of using either high 
speed computers or of finding an approximate solution; 
the latter course was chosen here. The problem is 
particularly suited to finding an approximate solution 
since certain relations (14), (15) have already been de- 
rived which the exact solution must obey; the degree 
to which they are satisfied by the approximate solution 
can then serve as a check on its accuracy. To review, 
relation (14) gives the melting rate which occurs when 
s = / (complete melting of the slab) while relation 
(15) gives the time after melting begins at which the 
slab is completely melted away. 

On examining Eqs. (17) and (18) we note that at 
z = 0 the function ¢ does not vary and hence 0¢ Or, -0 
= (. Since at the critical surface where melting actu- 
ally takes place O¢/Or is identically zero, an examin- 
ation of Eq. (17) suggests that we neglect this term in 
the equation during a first attempt at an approximate 
solution. If this is done 7 will no longer appear as an 
independent variable in the equation but will then only 
occur implicitly as a parameter. Consistent with this, 
we note further that on the boundaries of the region 
0<z< 1— Sthe function ¢ satisfies two time independ- 
ent conditions (18b, c). By neglecting O0¢/0r7 in Eq. 
(17) a linear second-order equation in z results which 
may be used to satisfy these conditions. The equation 


in question is 


vs oe (0/07) 79'(1 a; t*)¢ _ 

2? Oz T,/Q1 - S, t*) 
Dy AO io, 
Ty’(1 — S,t*) | TM(1— S, t*) 


Ce scl = 3, <«->0 (19) 


The solution [satisfying conditions (18b, c)] of this 
linear second-order differential equation may be ob- 
tained in general using the method of variation of 
parameters and is shown next. It will usually prove 
easier, however, to solve particular cases directly rather 
than to substitute into this general solution. 


o(z, 7) = Gle*— e*] —_ 


._ (% Fo(a, 7) e~*'* 
Az F 1 
f (Oe — ¥) da + 
*Fy ’ — Ara 
eu? {" dh alle (20) 
0 (re — d1) 
where 
S js? (0 ‘Or) To’ (1 —= =. f “ye 
A,pe=—--+ ——_———— 21 
, 2*\at na-sr (21) 
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HEAT CONDUCTION 


and 
G = l 7 
[1 — (Ne/Ales~G-S)} 
{ie Q, mes da — - e-MI=S) yg 
-S Fi(a, r)e~*** =a-S) 
J, . —s da + i, + (22) 


We shall denote by ¢, the first approximation given 
above by the solution of Eq. (19), to the function ¢. 
An improved solution may be obtained by forming 
0¢) /07r, now a known function of (z, 5S, S), and placing 
it in the inhomogeneous right-hand side of Eq. (17). 
In general, the 7th approximation to the function ¢ will 
be obtained by solving the equation 


0" ¢; 45 Oy; > (0/Or)70’/(1 — S, t*) os 

ds? | Og Pai=-sr 
0v.-1 To" (z, 7) STo'(z, t*) 
Or ia~- 2A” ht -ah 


Fide, ¢= 28... :< Ce 
satisfying conditions (18b, ¢c). This solution could 
again be written in general by the method of variation 
of parameters and would correspond in form to Eqs. 
(20)-(22). 

The process for approximating ¢ may be carried as 
far as is desirable. Note that the function ¢; deter- 
mined in this manner is still expressed in terms of 
S(r) and its derivatives and that one boundary condi- 
tion (18a) remains to be satisfied. When the desired 
approximation has been obtained the remaining bound- 
ary condition (18a) yields a nonlinear ordinary differ- 
ential equation for the determination of S(r). 

To summarize, instead of solving the original non- 
linear partial differential Eq. (17) our problem for de- 
termining the temperature in the solid and the amount 
of material melted at any time has been reduced to 
finding the solution of two ordinary differential equa- 
tions as follows: (1) the determination of ¢,(z, 7) in 
terms of S(r) and its derivatives to the desired degree 
of accuracy from a linear ordinary second-order differ- 
ential equation in z, and (2) the determination of S(r) 
from the nonlinear ordinary differential equation re- 
sulting from satisfying the third and last boundary 
condition (18a). It should be emphasized that in find- 
ing the solution of Eq. (23) the aid of high-speed com- 
puters is not needed. The method developed in this 
paper has the utility of readily providing accurate ap- 
proximate solutions using only a desk computer. 

We will use this method to find S(7) for two heat in- 
puts, the first constant and the second corresponding 
to that which might be encountered by a body descend- 
ing through the atmosphere. The constant heat input 
example will first be used to discuss the accuracy of the 
method by comparing several approximations with 
each other and with the exact relations (14) and (15). 


IN A MELTING SLAB 


te 
od 


Example I—Constant Heat Input 


The premelting (0 < ¢ < ¢*) temperature distribution 
in a slab (0 < x < /) under a constant heat input Q) at 
x = Oand insulated at x = /is'* 


Qol fxt | 3x? — 6xl + 21? 
k \P 6/? 


9 @ (—n?x*«t)/P? 
Yd “oe 


T(x, t) = 


nwx | 
cos F 
Ls 


9 9 


Tw 19 n* 


To obtain the postmelting solution we proceed in the 
manner shown in the previous section. However, for 
simplicity of illustration we shall assume that quasi- 
steady conditions have been achieved before the initi- 
ation of melting in which case the exponential terms 
may be neglected, though clearly this assumption is by 
no means essential to the use of the method. Expressed 
in the nondimensional variables of Eqs. (10) and (11), 
we then have at ¢ = ¢* 
To(z, t*)/T* = 1 + [(2? — 22)/2r] (24) 
Corresponding to Eq. (15) we find the nondimen- 
sional time after melting has begun at which the slab 
has been completely melted away to be 


7, = (1/3) + (n'/?/2)(r/M) (25) 


First Approximation 


We will now find the first approximation for the gen- 
eral case of arbitrary \/. We obtain from Eq. (24) 
To’ (1—S,t*) = S (0/dr)7%’(1 — S, t*) $ 


’ 


:" r ioe r 
To’ (z, t*) 2-1 To" (z, t*) l 
7 .* ig ae 
Using these relations in Eq. (16) our expression for the 
temperature then becomes 
O(z, r) = 1 + [(2? — 22)/2r] + (S/r)¢(z, r) (26) 
where the governing equation (19) for the first approxi- 
mation to the function ¢ is 


(02y,/dz?) + S(O¢,/dz) — (S 5S), = 
[((S — 1)/S] — (S/S)z (27) 


The solution of this equation for ¢)(z, 7) satisfying con- 
ditions (18b, c) is given by Eq. (20) or may be obtained 
directly as 





¢i(z, r) = A + Zz + Ce + CG™ (28) 
where 
Ayn2 = —(S/2) + 1 (S? ‘4) + (S/S)}12 
A = S-1+(1/S) 
(Az ,)Ae™—WU-S) 
= Th Cape RT] (28a) 
—A 

G= 


© [L = Qe/Ae™—MA“ 9] 


The equation resulting from substituting Eq. (28) 
into the remaining boundary condition (18a) for the 
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FIRST APPROXIMATION 





NON-DIMENSIONAL THICKNESS MELTED 


° 2 4 6 8 10 
a 
(Tio 
TIME SCALE 
Fic. 2. Nondimensional curves of thickness melted versus time 
for a constant heat input. Normalized scales. 


determination of (7) is 


\/2y 
stk 2 > 
QQ 2M 
(se— 1) (1—S) 
s? sa lt a9) 


1 — zt 
) {1 Zz (Az AY ated ied | 


We can already verify at this state that one of the de- 
rived, auxiliary relations, Eq. (14), which must be 
satisfied by the exact solution is obeyed by the approxi- 
mate one. Remembering that we are treating the 
case Q(r) = Qy, it is evident that at S = 1 we obtain 
from Eq. (29) the velocity of melting 


S|,a1 = 2M/r'r 


which is the melting velocity required in the exact solu- 
tion by relation (14). 

Eq. (29) is easily solved numerically by desk com- 
puter. The procedure used was the following: if S 
is known at some time 7;, the value of S(7;) can be 
obtained from Eq. (29) with Q(7;) = Qo; once the 
melting velocity S(r;) has been obtained the value of S 
at a later time 7;+; can be calculated using central dif- 
ferences—i.e., from 


S(rigt) = S(tin) + 2ArS(7;) (30) 


The process is then continued until S = 1—i.e., 
complete melting has occurred. 

This procedure is satisfactory everywhere except at 
7t = 0 where it fails for the following reason. We note 
first that setting S(7;) = 0 in Eq. (29) requires that 
S(7;) = 0 (for all L ¥ 0) and hence, if S(7;-1) also is 
zero, we find that S(r) = 0 for all 7 > 7;. Now since 
S(r) = 0 for r < 0 we have S(0) = 0 [ef. Eq. (4)]; 
hence S(7,) will also equal zero, and therefore the result 
S(r) = 0 for 7 > 0 would follow. We must therefore 
start the numerical solution of this equation in a dif- 
ferent way—-namely, by finding an analytic solution 
which is valid in the neighborhood of r = 0. In this 





region it can be expected!! that S < S$ < 1. Under 
these assumptions Eq. (29) reduces to 


Or) Qo =1+ (7/2 2)(r M)S iets {s $3 2 
which has as its solution for Q(r) = Qo 
S(r) = (2/m"!2)(M/r) [73/2/(3/2)3/2] (31) 


a function which satisfies the assumptions just made 
on its behavior and that of its derivative in the neigh- 
borhood of the origin. 

We can now see that the initial condition (18d) will 
be satisfied—.e., 

lim 7)/(1 — S, t*)¢,(z, 7) = lim S¢i(z, rT) = 0 (32) 

70 70 

Eq. (31) can be taken to be the solution for S(r) in 
the region 0 < +r < 7;. For +r > 7 the numerical 
procedure outlined above was employed. 

The solutions obtained in this way from Eq. (29) for 
S(r) are presented in Fig. 2 by curves of S(7r) versus 
t/(r,a for M/r = 0.75, 2, 4, and 50; (ra de- 
notes the approximate value of 7; (nondimensional 
time after melting at which the slab is completely 
melted away) obtained by solving Eq. (29). A com- 
parison is given in Fig. 3 between (7,)a and the exact 
value of 7, [Eq. (25)] as a function of M/r. The 
amount by which these differ is a measure of the accu- 
racy of the approximate solution. It may be seen to 
be everywhere less than 15 per cent; in the range of 
greatest physical interest, namely, that of large 1/ 
i.e., a small latent heat of melting—the error is much 


less. 


Second Approximation 

The second approximation ¢2(z, 7) is obtained in the 
manner previously indicated. First, 

Oe1/Or = A + [Ci + MCiz]e* + [Co + oCrz]e 
is obtained from Eq. (28) and substituted into Eq. (23). 
The resulting equation 
(0°¢2/02*) + S(Og2/dz) — (S/S) eo = 

((S — 1)/S] — (S/S)z + A + 
[Ci + MCiz]e*” + [Cy + hoCoz Je** (33) 
is then solved for ¢o(z, 7) yielding 


g(z,7) = A+2e+ Ge*+ GQed* + 
Dz” + Doz?e* 4 F,ze” 4. Foze™* (34) 





where 
FIRST APPROXIMATION 
130, 
120} 4 ae | 4 4 4 a 4 i 
(Tig | 
° 2 el 
110} + | + + + 4 ? + 4 
| 
100+ + + + + + + 4 : 


Fic. 3. The ratio of the approximate to the exact value of the 
nondimensional time after melting begins at which the slab is 
completely melted away versus the physical parameters occurring. 
Constant heat input. 
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HEAT CONDUCTION IN A MELTING SLAB 
Cy = —(ro/M) Coe -9 — [1 — S)/aa] {Di[2 + (1 — S)] + 
D.[2 + rl — Se} — (A/a) {All + A. — S)] + Fefl +A — Spe” | 
: — 1 = (1 — S) 
(;: = a=s), | 4 1D,[2 + (1 — S)] + D[2 + A(1 — S)]} + 35) 
2 it = t, Ayew—™ a |  s di iD, | 1 J [ lj (390 


| 
i, \Fi{1 + A(1 — S)] + Fell + Al — Se oma | 


Dy, = —\iCi/2(d2 — Ax), De = doC2/2(A2 — Ax) 
Fi = —{[C, (A2 — Aa)] — [ArCi/(A2 — Aa)? 

Fy, = [Ce (A2 — 1)] — [AeC2/(A2 — Aa)] 

A = (SS/S*) — 1 + (1/8) 


Corresponding to Eq. (29) for the first approximation, the equation for the determination of S(r) is obtained by 
substituting Eq. (34) into Eq. (18a), resulting in 
Q(r) _ va gil? ¢ a S[1 — e®—0-5) | 
Qo 2M i~- err 
{rod — (Ae/M)S(L — S)[Dif2 + (1 — S)} + Deof2 + ro(1 — Syfe™—A-Sy 
(Me/Ma)SUFif 1 + (1 — S)}+ F of 1 + As(l — Spe -9 1} + Ss. — S)[D\{2 +401 — S)} + 
D{2 + ro(1 — Spe 9-9] 4 STRf1 + a1 — S)} + Fell + r(1 — Spe™*-"-9] — SiR, + KF] (36) 


xX 


It is seen that at S = 1 the second approximation 
also results in the correct melting velocity - 


+ 







Slo, = (2/2'!?)(M/r) 


° 
@ 


given in Eq. (14) for the exact solution. 

} We shall solve Eq. (36) numerically in the same gen- 
eral manner as Eq. (29) was solved. The only differ- 
ence results from the fact that § now appears in the 
equation along with S and S which occurred in the first 
approximation equation. We proceed as follows: 
knowing the value of S and S at some time 7; the value 


1)—FIRST 
} APPROXIMATION 


MELTED 


2)- SECOND 


° 
on 


™ 
~~ *50 


THICKNESS 


° 
. 


of S(7;) was obtained by trial and error so as to satisfy 
Eq. (36) with Q(7r) = Qo. Once S(r;) had been ob- 
tained the value of S and S at some later time 7,41 were 
calculated using central differences from P 


02+ 


NON-DIMENSIONAL 


) oO! o2 03 o4 os ‘ 06 


S( 7441) =z S(7; 1) + 2Ar SMrit (37) TIME SCALE 
S(rigt) = S(rin) + 2Ar S(ri)§ Fic. 4. Nondimensional curves of thickness melted versus time 
for a constant heat input. 
and the process continued until S = 1; i.e., complete 
melting. 
To start the solution to Eq. (36) it is still necessary 
to find an analytic solution which is valid in the neigh- 
borhood of 7 = O since S(O) = ~. For S<S< 1 


of S versus 7 which are plotted in Fig. 4 for M/r = 
4.89 and 50 and compared with the corresponding first 
approximation curves. We note that the two approxi- 


about r = 0, Eq. (36) becomes mations for M/r = 50 are almost identical and that 
for M/r = 4.89 the approximations are nearly the same 

Q(r)/Qo = 1 + (4"/?/2)(r/M)S — over a major portion of the melting time. Further we 
(9 8), S St! stoves (S25 Ss){S S}} * + see that the error in (7,)a for M//r = 50 is negligible and 

(5/2S){ S/S} */? + (8/8S){ S/S}? that for M/r = 4.89 it has been reduced from 9 per 


which has as its solution for Q(r) = Qo. cent to 3 per cent from the first to the second approxi- 


mation. 
" \ ‘ 1/2 a 1/2 /J\3/2 3/2 22 *4° . ° 
S(r) = (2/a"?)(M/r) |(5/4)/(38/2)*!? |r (38) In addition, curves of S versus 7r/(7,)a are shown for 
a . : each approximation in Fig. 5 and Fig. 6. It can be 
and is of the same form as the starting solution (31) ; : , 
; : : ; seen that when normalized in this manner the two 
for the first approximation. ; ; 
: ; ee - aa approximations fall very close to each other over the 
We can again see that the initial condition (18d) is é : ial 
“eA ‘ total melting period. This suggests that a good ap- 
satisfied—1.e., : : : 
proximation to the exact solution for the amount of 
lim 7) (1 — S, t*)¢(z, 7) = lim S¢o(z, 7) = 0 (39) material melted at any time may be obtained from the 
vase — first approximation for S versus r by multiplying the 


Using Eq. (38) to start the solution we obtain curves ordinate by 7/(r,)a. 
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Fic. 5. Comparison of first and second approximations of the 
nondimensional thickness melted versus time for a constant heat 
input. Normalized scales. 
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The temperature in the slab at x = //2 and/ (z = 
(1/2) — Sand 1 — S) can be obtained from Eq. (26). 
For the particular case of M//r = 4.89, r = 0.75, both 


the first and second approximations are shown in Fig. 7 
as a function of 7/(7,)a. One sees that, as was the case 
for the thickness melted, there is little difference be- 
tween the first and second approximations for the 
temperature in the slab. 


Example II—Atmospheric Heating 


In this section we will calculate the first approxima- 
tion for the problem of a slab insulated at x = / but 
now subjected to a time varying heat input Q(¢) on the 
other face. The particular heat input used is meant 
to correspond to one that which might be incurred by 
a body descending through the atmosphere and is given 
in Fig. 8. The maximum rate of heat input was taken 


as 200 cal./cm.?-sec., (736 B.t.u./ft.?-sec.) with a total 
heat input of 2,244 cal./cm.*, (8,258 B.t.u./ft.2). The 
form of Q(t) was arrived at in the manner shown in 
reference 15. 

To obtain an analytic representation of the temper- 
ature in the slab a function of the form Q(t) = R(e* — 
1) was used to approximate the actual heat input Q(/) 
over the first portion of the descent time. Once melt- 
ing begins the actual heat input can be used since the 
procedure is then numerical. The agreement between 
this simple representation of the heat input Q(¢) and 
the actual function Q(f) is also shown in Fig. 8 for the 
numerical values given in reference 16 and with 


O = 11.776(e% 52" — 1) (40) 


Using the expression for Q(t) the temperature dis- 
tribution within the slab for 0 < ¢ < ¢* (i.e., before 
melting starts) is derived in reference 16. It is 


Rl J cosh y~@/2(1 — x/I) 
T(x, i) = <9 jal 1/2 ‘ 
= oF ry ? a 
] (x — 1)? 
_"- _ “re + 
oO y 4 
de ra Kt/l* nme) 
> nx? + nxt —_ 1§ (41) 
1,2 arty 
where + ie 


Consider a steel slab with the following physical prop- 


erties: 


1 = 2ecm., «x = 0.0986 cm.?/sec., 
p = 7.8 gm./cm.? 

L = 8cal./gm., 7* = 1,500°C., 
c = 0.13 cal./gm.-°C. 


(42) 


From this we may calculate that 

M = (7! 2)ter™ L) os 216, T= Ri™ 1Qo — 0.375 
where the reference value of the heat input was taken 
to be 200 cal. /em.?-sec. 
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Fic. 7. Comparison of first and second approximations for the 
nondimensional temperature versus time for a constant heat 
input. Normalized scales. 
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Fic. 8. Melting due to a heat input incurred by a body descend- 
ing through the atmosphere. 


It can be shown from Eq. (41) using these properties 


and the heat input given by Eq. (40) that at ¢ = 19 
sec. the temperature at x = 0 will be 1,500°C. and 
melting will begin—.e., 7)(0, 19) = 7)(0, t*) = 7*. 


To facilitate computations we express the nondimen- 
sional temperature at this time, 

A(z, 0) = To(x, t*)/T* 
in terms of a polynomial. The temperature 0(z, 0) is 
accurately given by 


To(x, t*)/T* = Gz, 0) = 
1 — 2.66672 + 3.32672? — 


2.1809z* + 0.639424 (43) 


where 6(z, 0) has the same value and slope at x = 0, / 
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Fic. 9. A nondimensional temperature distribution at the time 
melting begins in a slab undergoing atmospheric heating 


and the same value at x = //2as6@(z,0). A comparison 
of 6(z, 0) and 6(z, 0) in Fig. 9 shows them to be indis- 
tinguishable for the scale used. 
Using Eq. (43) we obtain from Eq. (16) 
O(z, 7) = 1 — 2.66672 + 3.32672? — 
2.18092? + 0.639424 + 
1.236S(1 — 0.91048 + 2.0685S7)¢:(2, 7) (44) 


The function ¢;(z, 7) is then the solution of Eq. (19) 
containing the inhomogeneous terms obtained from 
Eq. (43). Eq. (19) becomes 


(07¢,/0327) + Se) dz) — (S S)qei = 
a — bz + cz? — dz* (45) 


where 


1 — 1.8209S + 6.20398? 
4 ~*~ “1 — 0.91045 + 2.0688? 


(2.15758 — 5.383) 
a 


(5.29348 — 6.2039) 


S(1 — 0.9104S + 2.068S?)’ 


~ S(1 — 0.91048 + 2.06882)’ 


(5.3838 — 10.5869) 
S(1 — 0.91048 + 2.0688?) 
(2.068) 

S(1 — 0.9104S + 2.068S?) 


The solution of this equation satisfying conditions (18b, c) is 


gil, tT) = A,e" + Ase” + By + Biz 4- Boz? + B;2* (46) 
where Ape = —(S/2) + {(S? 4) + (S S)q} 1/2 
A; = — Bo _ Ao 


e~MO—S) 1) Boe" *) + 1 — B, — 2B(1 — S) — 3B;(1 — S)*) 


eS 3 1 — (A2/Are™*—O-S) 
and B; = dS/qS, B2 = —(cS/qS) + (3SB:/q; 
B, = (bS/qs) + (S/q)[(6B3/S) + Bz] 
Be = —(aS/qS) + (S/q)[(2B2/S) + Bil 


The third and last boundary condition for the determination of S(7) is given by Eq. (18a). 


this becomes 


In our present problem 
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Notice that at S = 1, the melting velocity is 


7 2 M Qr) Or) 
» = = = 65 
s=1 gl? ¢ Qo s=1 Qo gal 


which is once more the condition the exact solution 
must satisfy [cf. Eq. (14) ]. 


To start the numerical solution of Eq. (47) we must, 


again obtain an analytic solution which is valid in the 
neighborhood of zero. Assuming that S < S < 1 we 


obtain the equation 
Q(r)/Qo = 1 + (8/65) — 2.495}.S/S} 1/2 


for this region. Since Q(7)/Qo — 1 = 0 over the small 
region needed to start the solution we find that the re- 
maining terms are satisfied by 


S = 88.277 79/2, S§ = 132.4271 (48) 


functions which agree with the assumptions made on 
their behavior about 7 = 0. 

The solution for S(7) obtained from Eq. (47) and the 
heat input Q(r) it corresponds to is shown in Fig. 8. 
For the particular choice of variables used approxi- 
mately 3/10 of the slab was melted away. 
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Electrogasdynamic Motion of a Charged Body 
ina Plasma’ 
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(1) Introduction 


i has generally been concerned with 
the motion of bodies through compressible fluids. 
In recent years there has been renewed interest in a 
fluid which differs in so many of its properties from an 
ordinary gas that it has been called a ‘fourth state of 
matter.’ We refer, of course, to a plasma which is an 
electrically neutral gas containing ions and electrons. 
Plasmas now occur in many situations which are of aero- 
nautical or astronautical interest—e.g., the ionosphere, 
interplanetary space, jet exhausts, hypersonic shock 
layers, etc. The possibility of affecting the properties 
of plasma flows by means of external magnetic fields 
has stimulated an extensive literature in aeronautical 
magnetohydrodynamics. We wish now to consider 
the motion of a body through a plasma where electro- 
static phenomena are significant. 

The motion of body in a plasma will in general induce 
a negative charge on the body due to the fact that the 
relative flux of the electrons to the body is greater than 
that of the ions. This charge will then affect the flow 
through the Coulomb interaction with the plasma. 
A charge may also be induced on a body by entirely 
different or even artificial mechanisms. The motion 
of charged bodies in a plasma is what we have called 
“electrohydrodynamics,”’ or since a plasma is generally 
a compressible fluid. “‘electrogasdynamics.”’ 

The effect of an electric charge on the motion of a 
body through a plasma was discussed by Jastrow and 
Pearse! in connection with satellite motion through the 
ionosphere. They considered the effect of charged 
particle interaction with the static potential field of the 
charged body, assuming that the electrostatic field is 
undisturbed by the flow. However, a plasma flow 
interacts with an electrostatic field on the body, and 
the dynamic modifications of both the field and the 
This was done by Kraus 
considerable 


flow must be considered. 
and Watson.? The modifications are 
when the flow velocity is greater than the velocity of 
sound in the (ion) fluid, as is the situation for satellite 


or meteorite motion. 
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It is our purpose in this paper to review briefly the 
basis for treating a plasma as a fluid to which hydro- 
dynamic considerations may be applied. We will then 
describe briefly the basic equations and illustrate the 
theory by treating the very simple case of the super- 
sonic flow of a plasma past a two-dimensional slender 
dielectric body. 


(2) The Plasma as a Fluid 


A plasma differs from a gas of neutral particles in 
certain important respects. In a neutral gas, local 
properties are determined by collisions between par- 


ticles. The free path for such collisions is given by 


l = 1/mnq 


where 7 is the number density of particles and 7q? is 
the effective cross section of the particle where 4g, 
for example, will simply be the radius of neutral par- 
ticles when considered as rigid elastic spheres. In 
general, g is the distance of closest approach under the 
influence of the force exerted by one particle on an- 
other, the force here being so short in range that only 
two particles are involved in any encounter. 

In a plasma, however, the influence of the inverse 
square Coulomb force between charged particles is 
effective at considerable distances; thus, the motion of 
a given particle influences the motion of many par- 
ticles. Alternatively, each charged particle moves 
under the influence of all the charged particles in the 
With respect to a particular test particle, 
First, 


plasma. 
the plasma may be separated into two regions. 
there is a region in the neighborhood of the test par- 
ticle in which the effect of the Coulomb force is similar 
to that of short range forces, the interactions here 
being of a random thermal character. Second, there 
is the region exterior to this neighborhood, in which 
the individuality of the particle charges need not be 
considered and the overall effects of the Coulomb forces, 
due to these particles, can be considered as a body 
force on the test particle. The size of the region in 
which the random thermal motion can occur may be 
characterized by the Debye radius, 
R,* = kT /4xn,e* 

where 7 is the plasma temperature in degrees Kelvin, 
n, is the number density of electrons, e is the electron 


charge, and k is the Boltzmann constant. The Debye 
distance may also be characterized as the distance in 
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which electric charges can be separated by the thermal 
motion in the plasma. Since the test particle in a 
plasma is affected by all the particles in the Debye 
mq”, may be de- 
These are 


region, an “effective cross section,” 
fined, as well as an “effective free path.”’ 
similar to the quantities defined above for the inter- 
action of two neutral particles, except for a factor 
which takes into account the effects of all particles in 
the region.* + We may then define a distribution 
function which satisfies Boltzmann’s equation, and the 
collision integral now refers to interactions in the modi- 
fied sense just described. By taking moments of the 
Boltzmann’s equation, we arrive at the usual macro- 
scopic equations of motion. These equations can be 
reduced to the ordinary equations of gasdynamics 
when a scalar pressure can be introduced. This can 
be done if the ‘‘effective free path” is less than a char- 
acteristic length of interest in the problem. We note 
now that this is possible in the case of a plasma when 
the characteristic length is much less than with neu- 
tral particles. The effective cross section for charged 
particles may be several orders of magnitude greater 
than for neutral particles in a gas at the same density 
and temperature.* 

For the ionosphere and in interplanetary space, the 
ionization density is very low and the “effective free 
paths” are very long. It may not, therefore, be pos- 
sible to reduce the pressure tensor to a scalar. If, 
however, there is a charged body in this thin plasma, 
there will be a macroscopic electric field. There will 
now be a collective effect similar to that produced by a 
point charge. This effect now appears in the Boltz- 
mann equation as a macroscopic body force. By 
introducing the directed particle velocity resulting from 
the electric force and an approximate linearization 
based on the difference between the directed and ran- 
dom velocity, the equations describing the flow can be 
reduced to a form (at least when dissipative effects are 
neglected) which is analogous to those for the linearized 
continuum flow.” The physical basis for treating a 
rarefied gas with an impressed electric field is similar 
to that encountered in the magnetohydrodynamics of 
rarefied high temperature plasmas, but there is con- 
siderable simplification in the former because the 
complicated electromagnetic effects are absent. 


(3) Basic Equations 


To exhibit some of the phenomena of electrogas- 
dynamics, we consider the flow equations which, in 
view of the above discussion for the high plasma density 
case, may be written in the steady case as 


MnN-VV + Vp; — neE = 0 (1) 
V- (nV) =0 (2) 
pm,’ = const. (3) 


where ./ is the ion mass, m; the number density of ions, 


* For a more detailed discussion, see reference 3, Chap. 5. 
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p; the ion fluid pressure, y the adiabatic constant, 
V the velocity vector of a fluid particle, and E the elec- 
A similar set of equations governs 
Here, it has been assumed 


trostatic force field. 
the motion of electrons. 
that the plasma is composed of electrons and singly 
ionized atoms of one chemical species. The electric 
field is determined by Poisson’s equation: 


V-E = 4ne(n; — n,) (4) 


If the electrons are in local thermodynamic equilibrium, 
it is readily shown from the equation of state and the 
electron equation of motion with the inertia term 
neglected that® 


Ne = N, exp (e&/RT) (5) 


where V@ = —E defines the electrostatic potential and 
n, is the equilibrium ionization density. Furthermore, 
the electron pressure », may be found from the electron 
equation of state, p, = n,k7T. or using Eq. (5), we ob- 
tain 


pe = n kT exp (e&/kT) (9a) 


(Note that the equilibrium electron temperature 7 will 
be constant throughout the flow.) 

If we now define the usual small disturbance vari- 
ables by the formulas 


Po + Pir 
V =_ Vo + v 


N= +n, t= 
nN, = n,(1 + e®/kT), 


(where ,, ~,, and Vp refer to the undisturbed free 
stream plasma, and n’, p,’, and v are small perturba- 
tions), and we introduce a velocity potential P by 


VP =v (6) 
we obtain the coupled linear equations 
—V,*P,, + a?V?P = (eV,/M)®, (7) 
[V? — (A*/Rp*a;”)]® = (wp,?/ai?)(V.M/e)P; (8) 
Here, we have defined the sound speed in the ion fluid, 
a;, by 
a, = yp./nM (9) 


x is the coordinate in the direction of the free stream, 
Rp is the Debye distance defined previously, and 
A? = a;? + w,/*Rp’ (10) 
where the quantity 
Wp; = 4mre*n,/M 

will be recognized as the ion plasma frequency. 

Let us next introduce the following scale transforma- 
tion in Eqs. (7) and (8). 


/ / 
E=x/c, n= Y/C, 


f= s3/c 
P(x, y, 2) = cV,p(&, 7. £) ®(x, y, 2) = 
(MV,?/e)¢ (& 1, 6) (11) 


where c is a characteristic length of the problem, and 
(x, y, s) are Cartesian coordinates. 
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ELECTROGASDYNAMIC 


The basic Eqs. (7) and (8) now become 
aaa Vo? Pe + a;,V*p aad V7: 
[V2 — (A%c?/Rp*a;”) |]o = c*(wp;?/a;7) p; 


(12) 
(13) 


Eqs. (12) and (13) may be decoupled to obtain, for 
example, for p, the fourth-order equation 


[V2 — (A°c?/wy,?Rp*)| X 


(—V2be + 4°V*p) — (CV 2ap?/a) pe = 0 (14) 


A similar fourth-order equation may be obtained for ¢. 
The equation of the characteristics of Eq. (14) is a 
fourth-order polynomial. If the free stream velocity 
V, is less than the ion sound speed, all the roots are 
imaginary, and the equation is completely elliptic. If 
the V, is greater than the ion sound speed, we have 
two real and two imaginary roots. That is, the equa- 
tion is twofold elliptic and twofold hyperbolic. The 
latter case is of particular interest, since here one ex- 
pects many of the effects of supersonic aerodynamics. 
For an electrostatic point source, a solution of this 
equation has been obtained.” The solution shows that 
for the velocity of the source less than the ion sound 
speed, no significant effects are obtained, except that 
the electrostatic equipotential lines are distorted from 
their usual spherical shape to ellipsoids. For a source 
velocity greater than ion sound speed, these equipo- 
tentials extend to infinity along the characteristics. 
A “bunching up” of equipotentials occurs near the 
Mach cone; upstream, the potential decays exponen- 
tially. Downstream of the Mach cone, the equipoten- 
tials decay algebraically, and there is an important 
transition region where the ionization density shows a 
In addition, there is a supersonic 
drag phenomenon. It should be noted that essentially 
the same solution was also obtained for the case of very 


marked increase. 


low density plasmas. 
The Debye radius is usually very small in comparison 
with other lengths of interest in the problem, so that 


V°o/¢ << A®c?/a?Rp? (15) 
Eq. (13) now reduces to 
@ = —(1/V,?)M,’Ro*wp,"p; (16) 


where M, = V,/A is the “plasma’’ Mach Number. 
Eq. (12) with (16) now becomes 


—M,2pe + V2p = 0 (17) 


The reduction of the order of the equation by this 
singular perturbation has important consequences. 
The second-order equation will be valid throughout 
the flow, with the possible exception of certain regions 
such as the neighborhood of the body where boundary 
conditions are required. In these regions where the 
inequality, Eq. (15), does not hold, there appears a 
type of boundary layer phenomenon. It will also be 
noted that the Mach angle is altered from the exact 
case with the acoustic velocity increased to 


A = (a? + wp?ZRp’)” 


MOTION OF 


bo 
we) 
— 
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from the exact value a;. The resulting steepening of 
the Mach waves may be interpreted in the following 
way. The twofold ellipticity of the fourth-order 
equations implies that perturbations of the flow at a 
given point will propagate upstream to infinity, the 
intensity of the perturbation attenuating as it propa- 
gates upstream. In the reduced system, this upstream 
propagation will, however, be confined between the 
region formed by the two Mach cones resulting from 
the two acoustic velocities. The elliptic character of 
the flow in this region assumes a modified hyperbolic 
character. In the exact case, there is, of course, addi- 
tionally a weak effect upstream of the steeper Mach 
cone. We have essentially then an ‘electrostatic 
Mach cone.” 

To define a unique solution to the problem, we must 
supplement the above differential equations by bound- 
ary and Cauchy type conditions. The number and 
type of conditions required to yield a_ well-posed 
problem follows from the purely fluid or electrostatic 
problems. These include the following. At the sur- 
face of the body, the direction of the velocity vector 
must be parallel with the slope of the body. For the 
electrostatic potential ¢, we have different conditions 
at the body depending upon the nature of the body 
We shall consider it composed of either a 
For the con- 


material. 
conductor or an isotropic dielectric. 
ductor, the potential must be a constant at the body 
surface; this constant is related to the total charge on 
the conductor, since the normal derivative of the po- 
tential at the body surface is essentially the local sur- 
face charge density. For the dielectric, the potential 
or its normal derivative may be arbitrarily prescribed. 
The surface charge accumulation may be determined 
in this case from the plasma flow problem. together 
with the purely electrostatic problem in the interior 
of the dielectric body. This will be illustrated in a 
subsequent example. 

Far from the body, the function p and the potential 
@ must either vanish identically or be composed of 
“outgoing waves,” as in the usual sense of supersonic 
flows. 

For the decoupled fourth-order equations such as 
Eq. (14), clearly all of the above conditions, expressed, 
of course, in terms of the decoupled dependent vari- 
able, must be used. 

Lastly, we obtain from the linearized equations of 
motion of the ions the ‘‘Bernoulli equation” 

pi’ = —n,MV,? (ph; + ®) (18a) 

For the electrons, we obtain from Eq. (5a) the ex- 
pression 

pb.’ = n.MV,2 (18b) 


where p,’ = p, — pp. 


(4) Plasma Flow Over a Thin Dielectric Body 


As an example, we shall consider the two-dimensional 
flow over a thin profile at zero angle of attack. To 
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Fic. 1. Electrogasdynamic motion over a dielectric body. 


simplify the mathematics, we shall choose an example 
for which the approximation, as given in Eq. (15), 
is valid throughout the domain of the problem. This 
will mean that “boundary-layer effects’’ must be mini- 
mized. The body selected must thus have cusped 
leading and trailing edges to avoid the ‘boundary 
layers’ along the nose and tail Mach waves. Addi- 
tionally, no condition on ¢ will be directly imposed 
at the body. Instead, ¢ will be so selected that no 
boundary layer arises along the body. 

The symmetric profile to be considered is given in 
the upper half plane by the equation 


n = (7,/4)(1 + cos 2rt) — (1/2) < € < (1/2) (19) 


where 7, is the thickness ratio. The profile is quali- 
tatively sketched in Fig. 1 where also the orientation of 
the coordinate system is given. Our basic flow equa- 
tion can now be written as 


where \* = J/,? — 1. This decoupled equation for p, 
together with the conditions on p as given previously, 
will be sufficient to determine a unique solution. This 
will be analogous to the Ackeret solution. The result- 
ing p will be inserted in Eq. (16) to obtain ¢. We ob- 
tain to this order of approximation 


p = —(1,/4d) cos 27 (E — An), — 3 < E-— Mn <3 
& = }[r(d2 + 1)]/2AV,2! Rp’wp,27, X } 
sin 27(~ — An) J 
p= 0, §-dAy< —3, E—A>}3 | 
$= 0 
(21) 


Note that p and ¢ are symmetric with respect to 7. 

From Eq. (21), the drag may now be computed. It 
will be composed of the fluid pressure drag (wave 
drag) Dy and the electrostatic body force drag, Dz. 
The calculation of the wave drag may be calculated by 
integrating the pressure, as given by Eq. (18), over 
the contour. The result is 


Dy = (m?/4)(m0M V,?7,2c/d) (22 


or 
Cop = (4?/2)(t62/d) (23) 


Here and in the following, we have assumed that the 
electron and ion temperatures are equal. 
The electrostatic drag Dz is given by 


+(1/2 
De = 2(M/e) f  @(€)QG/d#)dE (24) 
— (1/2) 


where w(£) is the induced surface charge density given 
by® 


€,(O@ On) oad €p(Odz, On) = — (€¢ MV,?)w (€) (25) 


Here, ¢, and eg are the dielectric constants for the fluid 
and the body, respectively; ¢, is the potential in the 
interior of the body, and the normal derivative to the 
body surface is to be taken as positive in an outward 
direction from the body. 

The function ¢, is now determined from the purely 
electrostatic problem in the interior of the dielectric 
body. In the absence of volume charges, ¢, is defined 
by the Dirichlet problem where the boundary values 
are given by Eq. (21) with 7 set equal to zero. 

Assuming a slender body, one may solve this prob- 
lem approximately with the result 


dn = —{[r(d? + 1)]/2AV,2} r.Ro'wp?? X 
sin 27é cosh 27 y= (26) 


This expression, together with Eq. (21) for ¢, may 
now be used in Eq. (25) to compute the surface charge 
density. Integrating the charge density over the 
body now gives the total charge on the body. In the 
present case, the net charge is zero. The electrostatic 
drag is next found from Eq. (24) and is given by 


Dg = ('*/d)(M/e)*(€,/c) Vo'7,?[1/(y + 1)7] = (27) 
The total drag coefficient is now given essentially by 
the sum of Eqs. (27) and (23), that is 
Cp = Cor (1 + 7) (28) 
where 


r = Dzg/Dp = (160°/d)M,*(Rd/c)*[1/(y + 1)] 


Here we have set «, = 1, and we have assumed equal 
electron and ion temperatures, so that 


9 9 9 
y= Ywp;?Rp* 


Since for a dense plasma, Rp/c < 1, we see therefore 
that the electrostatic body force is negligible in com- 
parison to the wave drag. 

On the other hand, it is seen that the wave drag in 
the plasma is significantly different as compared to the 
neutral gas case at the same thermodynamic state. 
The ratio R of the plasma and neutral gas wave drags 


being given by 
R = VM,? — 1/VM,? — 1 


where J/,, the “plasma Mach Number”’ (with A? = 
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a;° + w,?Rp? used as the square of the acoustic speed), 
is ¥ ‘y (y + 1) times the Mach Number J/, based on 


the ion thermal speed. 


(5) Final Remarks 


In the present paper, a simple example has been 
examined which indicated some of the aspects of the 
coupling between the fluid motion and the electro- 
static properties of a plasma. The simplication re- 
quired that “boundary layer” effects be absent. Here, 
“boundary layers’ are defined as narrow regions in 
which the electrostatic potential and the flow proper- 
ties change rapidly. In a more general case, important 
modification must therefore arise near the body. The 
precise analysis of these effects would require the use 
of the complete fourth-order system of differential 
equations. As noted previously, these equations are 
twofold hyperbolic and twofold elliptic when the motion 
of the body is supersonic with respect to the ion thermal 
velocity. This property manifests itself, for example, 
by an upstream effect ahead of the leading edge Mach 
wave which damps exponentially. Electrical effects, 
instead of influencing the entire medium, are confined 
essentially to a region downstream of an electrostatic 


Mach envelope. Along these Mach envelopes, charge 


separation takes place which provides an interesting 
modification of the flow that may be observed experi- 
mentally. 

The dense plasma case treated in the present paper 
is, of course, not directly applicable to the motion of a 
charged body in the rarefied ionosphere. In the latter 
case, ‘fluid type’ equations also arise when the molecu- 
lar collisions due to short range forces are neglected, 
and a treatment similar to the dense case may be car- 


ried out. 
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A Study of Vortex Cancellation 


(Continued from page 196) 


Moreover, the value of s/a is given by elliptic 
Combining this result 
0.62 for vortex can- 
Thus, the analysis indi- 


(6.11). 
wing theory? as s/a = 6.45. 
with Eq. (6.11), one finds Mmar 
cellation in a circular duct. 
cates that for the geometry considered only 62 per 
cent of the pressure drop.across the first airfoil can be 
recovered even if there are no skin-friction and separa- 
tion losses in the system. This value can be compared 
to the experimental recovery efficiencies of 50 per cent 
that have been realized. 

In concluding the discussion, mention will be made 
of two interesting side calculations that can be per- 
formed using the results of this study. The first calcu- 
lation involves the ratio Ap,/Ap;. If one uses the ex- 
perimental value 7 = 0.5 and solves Eq. (6.3) for the 
above ratio, there results Ap,/Ap; = 3. This implies 
that the induced drag comprises 3/4 of the total drag 
The second calculation involves the 
Since, by Eq. (6.10), Api)min = 


on the airfoil. 
ratio Api) min Api. 


Jow, one can evaluate AP;)min experimentally from tuft 
grid pictures. Moreover, Ap; can be measured in the 
manner previously described in the discussion concern- 
ing the experimental determination of 7. It is found 
that for the tests performed AP;)min/Ap; ~ 0.8. This 
implies that 4 5 of the irreversible pressure drop across 
the airfoil comes from establishing the vortex field and 
1/5 comes from skin-friction and separation losses. 
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On Maximum Pressure Rise With Turbulent 
Boundary Layer 


A. Mager and A. W. Goldstein 

Research a, ASTRO, The Marquardt Corporation, 

Van Nuys, Calif., and Aeronautical Research Engineer, NASA, 
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November 5, 1959 

a B. S. Stratford published’? the results of an 
analytical and experimental investigation concerning the 

maximum attainable pressure rise in the presence of incom- 

pressible turbulent boundary layer. The examination of his 

experimental data reveals that it is also in essential agree- 

ment with some of our work,’ in which the allowable turbulent 

boundary-layer pressure rise was used as a criterion for the pre- 

scription of cascade blade velocity distribution. A comparison 

between the analytical results of reference 3 (computed, how- 

ever, for the diffussion coefficient D7 = —0.017) and the Strat- 

ford experimental data is shown in Figs. 1, 2, and 3. 

As may be seen from Fig. 1, the agreement of the analytically 
predicted era velocity distribution is not too good initially, 
but, for (U/U;)? < 0.65, the data a1e very close to the results of 
our calculation aical 
from the fact that the turbulent boundary layer can withstand 
some suddenly applied pressure rise (such as produced by shocks, 
for example) without separation, while our calculation scheme 
In this respect, the analyti- 


This initial discrepancy stems probably 


did not allow for such happenings. 
cal method suggested by Stratford is somewhat better in this 
region (see Fig. 4 of reference 2). It should be noted, how- 
ever, that, from the practical standpoint, this difference in the 
velocity distribution in the initial region is rather unimportant, 
since it occurs for only very small surface length (1.0 < x/x» < 
1.15) and, thus, should not affect substantially the aerodynamic 
design. 

Turning now to Fig. 2, the predicted behavior of the momen- 
tum thickness (0/0) as a function of the pressure coefficient Cp 
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[= 1 — (U/U,)*] shows good agreement initially but tends to 
overestimate somewhat the boundary-layer momentum defi- 
ciency for very high pressure coefficients. As will be mentioned 
later, this probably is the effect of an outmoded friction law 
which our early method has employed. On the other hand, the 
development of the form factor H, shown in Fig. 3, while ini- 
tially slower by calculations, rises more rapidly later on to reach 
ultimately nearly the same final values as the experimental 
data indicates (H = 2.6). This behavior of the form factor H 
reflects, of course, the somewhat steeper initial experimental 
pressure gradient, which may be expected to affect primarily the 
shape of the boundary layer velocity profile. Even the 
predicted values of H are never too far off from those obtained 
experimentally and are certainly considerably better than those 
resulting from Stratford’s analysis (see Fig. 6 of reference 2) 
which underestimates the final value of H by a rather substantial 
amount (H = 2.0). 

It may be worthwhile to point out that this reasonably good 
agreement between the experiments of Stratford and the analysis 
of reference 3 is even more remarkable when one considers the 
fact that the analysis of reference 3 was based on a now out- 
dated correlation of Garner,‘ which was obtained from only 
limited amount of data, all measured with pressure distributions 
radically different from those suggested here. Furthermore, 
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Garner’s correlation resulted in a single analytical statement valid 
for all unfavorable pressure gradients and employed a zero pres- 
sure gradient friction law, both of which are thought today to 
be totally unacceptable for turbulent boundary-layer analysis. 
One cannot help but wonder, therefore, whether the agreement 
shown here indicates that the modern turbulent boundary-layer 
solution refinements are only a small perturbation on the essen- 
tial features already conceived by Garner and other early in- 


vestigators. 
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Pitot-Probe Measurement of Local Skin- 
Friction Reduction Due to Transpiration 
in a Low-Speed, Turbulent Boundary Layer 


Bernard M. Leadon and E. Roy Bartle 
Senior Staff Scientist and Staff Scientist, Respectively, 
Convair Scientific Research Laboratory, San Diego, Calif. 


October 13, 1959 


. Is perhaps unnecessary to mention that a fixed pitot probe 


set at about 
yor = yu Cr 2/v = 100 


is capable of providing the value of the local skin-friction co- 
efficient for zero blowing, C;,, by means of the usual law of the 
wall (see Table IV of reference 1). This well-known possibility 
has not been particularly recommended because, for most solid- 
wall situations, a simple round probe set firmly against the sur- 
face may be used equally as well.2 *7 However, the addition 
of a rigid foot to the probe will probably often be a permissable 
complication, and one may then freely use the law of the wall 
for Cy, and seek to extend the idea of such measurements to the 
transpiration case. 

Recently, one of us observed that the logarithmic portion of 
the low-speed velocity profile when expressed in certain coordi- 
nates has a slope which varies but little, if at all, with trans- 
piration rate.‘ As long as this is true, an exceedingly simple 
measurement of the local turbulent skin friction, with blowing 
or without, can be made. 

The logarithmic portion of the mean velocity profile for a low- 
speed, turbulent boundary layer on a porous flat plate may be 
represented by the equation 
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a ™% TN 4017 

a ¢ | u/u, AT y=.040" 

PAN 2 «/ u,=50 FPS ‘a 
w| v,/u,= 1,2,3,5 x10 
NJ 
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' —— 
CG % 56 LOG,(C,/C, ) = (u-u,)/uVG 72 
“4, oN —— REF. 4 
b sr) 
2 ig 
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1+ 2 
Cyu, 
Fic. 1. 
u+/C = A log(yt/C) + B (1) 


where 
a+ = u/uV C;/2, yt = yur C;/2 /v 


subscript 1 indicates the free-stream value, C is a function of the 
blowing rate having the value unity for zero blowing, and A and 
B are constants which may be taken as 5.6 and 4.9, respectively. 
When the available experimental data® are plotted in the coordi- 
nates u/u,; vs. log (y u/v), it is found that the blowing rate does 
not noticeably alter the slope of the logarithmic portion. Hence, 


1/C? = C;/Cyo (2) 
and, for any fixed value of y in the logarithmic region, 


A logi(Cy/Cyo) = (4 — to) /t1V Cro 2 (3) 
where subscript 0 indicates zero blowing values. The tendency 
of the logarithmic region is to expand with blowing, so that, for 
a fixed probe centered in the zero-blowing logarithmic region, 
Eq. (1) will probably remain valid up to quite high blowing rates. 

Alternatively, the probe may be moved so as to hold u/ con- 
stant, in which case 


Cr/Cyo = Yo/¥ (4) 


but there is some danger that, in moving the probe to larger dis- 
tances from the wall as blowing is increased, it may leave the 
logarithmic region and thus invalidate Eq. (1). 

The data of Mickley and Davis for the probe positioned at y = 
0.040 in. and for u, = 50 ft./per/sec. were used in Eq. (3) with 
A = 5.6 to compute the ratio C;/Cyo. These results are plotted 
in Fig. 1, together with an empirical curve which well represents 
the results of a similarity analysis of the complete velocity pro- 
files.4 The blowing parameter log (1 + 2v,,/Cysu,) most smoothly 
correlates such data, after which the results may readily be con- 
verted to the more practical parameter 27y/Cy,t1. 

It should be noted that some of the data scatter of Fig. 1 is 
surely due to the fact that the probe was repositioned for each 
reading, whereas, in an experiment planned for the use of Eq. 
(3), the probe would remain fixed. The contents of this note 
have been included orally with a related presentation.® 
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Limit Analysis of Multicell Boxes Submitted 
to Pure Torque 


Fernando Venancio Filho 
Instituto Tecnolégico de Aeronautica, Sao Paulo, Brazil 
October 9, 1959 


a an ideal elasto-plastic material, expressions for the 
collapse torque of two types of two- and three-cell boxes 
are derived. Results show that the collapse torque is, at most, 
of the order of 10 per cent higher than the torque at which yield 
begins. 

Consider the box of Fig. l(a), submitted to pure torque. The 
three following assumptions are made: (a) the thickness ¢ of the 
walls is constant; (b) the material has an ideal elasto-plastic 
diagram in shear, Fig. 2; and (c) the walls yield before buckling 
can occur. A purely elastic analysis furnishes the following 
values for the shear flows g; and go, taking into account that 
A = b/a: 
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Fic. 2. Ideal elasto-plastic diagram in shear. 
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2 + 4.54r T 
= SS a eee (1) 
1.57 + 5.56 + 13.062 a? 
1 + 6.53 7 
qe = (2) 


1.57 + 5.56 + 13.06d? a? 


a 


In comparing Eqs. (1) and (2), it is seen that g: > gz for 0 < 
\ < 0.5, and q2> gq: for \ > 0.5. For the first of the above con- 
ditions, yield will start in wall 123 when 


qi = de = Set (3) 
where 3, is the shear yield stress of the material, Fig. 2. Sub- 


stituting Eq. (1) in Eq. (3), one finds 7,, the torque which first 
causes yielding: 
ee 1.57 + 5.56 + 13.06? 
l,= = gea? (4) 
2+ 4.54r 
Under this torque shear flow in cell 2 is, from Eqs. (2) and (4), 


| 


dre = [(1 + 6.58d)/(2 + 4.54A)] ¢- (5) 


For torques TJ greater than 7,, the difference T — 7, will be 
resisted only by cell 2 until, for a torque 7;, the outer walls from 
this cell will yield, causing collapse of the whole box. Then 
7T..can be calculated by 

[( T. — T.) 2ab] + dre = Ye (6) 

Inserting Eqs. (4) and (5) in Eq. (6) and noting that \ = b/a, 


T. = (0,78 + 2dr)qea? (7) 
Finally, dividing Eq. (7) by Eq. (4), 
I, 1.6 + 7.6 + 9.12 


— s eae (8) 
i 1.6 + 5.6 + 13.1d? 


For the second condition, yielding begins in walls 34, 45, 


and 51. The counterparts of Eqs. (3)—(8), respectively, are 
qz = Ye 
T, = 1.57 + —_ > 13.062 = 
1 + 6.53 
die = [(2 + 4.54A)/(1 + 6.53A)] Ge 
(Te — T.)/0.25ra?] + gic = Ge 
T. = (0,78 + 2d)q.a? 
and 


a 0.8 + 7.1\ + 13.1° 


-- (9) 


ie 1.6 + 5.6 + 13.1? 
The results obtained from Eqs. (8) and (9) are presented in 
Fig. 3. 
For the doubly symmetrical box of Fig. 1(b), the purely elastic 
analysis gives, taking into account that by symmetry the shear 


separate boxes 123 and 
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flows in walls 123 and 1’2’3 are equal and noting that A b/a, 
2 + 4.54 7 
“1 = (10 


3.14 + 11.154 + 13.06? a? 

2 + 6.53 7 

- (11) 
3.14 + 11.154 + 13.06? a? 


Eqs. (10) and (11) show that go> q:forX > 0. Thus, yielding 


first starts at walls 11’ and 33’ at a torque given by 


a 3.14 + 11.15 + 13.062 

T. = gea* (12) 
2 + 6.53 

For torques greater than 7, the box will behave like two 
’2’3’ until it collapses at a torque Ne 


given by 


_ 3.14 + 14.83 + 13.062 
I. = , gea* (13) 
Dividing Eq. (13) by Eq. (12), 
T. 3.1 + 14.3A + 13.1d2 
i 3.1 + 11.24 + 13.1,? 


(14) 


The values from Eq. (14) are plotted in Fig. 3. 


Inertia Effects on Some Static Aeroelastic 
Problems 


D. J. Johns 

Lecturer in Aeroelasticity, College of Aeronautics, 
Cranfield, England 

October 20, 1959 


SYMBOLS 


a = normal acceleration 

B = moment of inertia about center of gravity 

c = wing chord 

d = longitudinal acceleration 

D = total drag 

e = distance of body centers of pressure from datum (Fig. | 
E = Young's modulus 


fraction of chord of wing center of pressure forward of elastic 


axis 
g = distance of body centers of gravity from datum (Fig. 1) 
G = shear modulus 
h = fraction of chord of inertia axis forward of wing elastic axis 
I = second moment of area 
J = polar moment of area 
k = stability parameter 
Ke, K@ = joint stiffness, critical value of Kg 
l = wing semispan, body length 
L = lift/unit incidence 
m = mass/unit length of body 
M = total mass 
qd pitching acceleration 
s = rigid static margin 
T thrust 
x = streamwise coordinate (Fig. 1) 
& normal offset distance of deformed body center of gravity from 
body datum 
a = rigid incidence 
6 = twist deformation 
Subscripts 
Ww = wing 
7 = concentrated mass 
1,2,c¢ = body (1), body (2), rigid combination, respectively 
A = aerodynamic 
I = inertial 


DISCUSSION 
TATIC aeroelastic problems are defined as those in which 
frequency-dependent air and inertia forces are not con- 
sidered, but inertia forces due to overall body motions may be 
included. It is with these latter forces that this note is con- 
cerned. 
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Fic. 1. Schematic of rigid combination. 


For unswept wings where the section flexural axis and 
center of gravity are nearly coincident, it has been the usual 
practice to ignore the lateral and longitudinal inertia forces in 
the analysis of wing torsional divergence. While this has been 
justified for the lateral forces—i.e., those normal to the wing 
plane, the neglect of the longitudinal forces has not been justi- 
fied. The main phenomenon that can occur is flexure-torsion 
instability, in which loads applied in the plane of a slender beam 
can cause deformations out of that plane if the flexural stiffness 
in the normal direction is low. One deduces from reference 1 
that for a thin, unswept, uniform cantilever beam—e.g., a wing, 
with the flexural axis and center of gravity coincident, the inertia 
instability parameter, for a uniformly distributed and a concen- 
trated tip inertia force, is 

k, = (dl2/\/GJEI)\ Mw + 3.225My) = 12.9 (1) 

For wings with, say, the center of gravity well forward of 
the flexural axis, both lateral and longitudinal inertia forces 
should be more significant in the torsional sense, being stabi- 
lizing and destabilizing, respectively. The aerodynamic tor- 
sional-divergence parameter for a uniform unswept wing is 
given by? 

ka = Lw(fel/GJ) = 2.47 (2) 
The corresponding criterion when lateral and longitudinal inertia 
effects due to uniformly distributed and concentrated tip masses 
are considered is easily shown to be 
kar = (cl/GJ){Lwf — ((2Lw/M.) — d] X 

(Mwhw + 2.47Mrhr)} = 2.47 (3) 

The inertia effects on most practical aircraft designs should be 
small, but, for highly swept wings, lateral effects might be signi- 
ficant. 

The problem of aerodynamic chordwise divergence in super- 
sonic flow is well known and need not be discussed further here. 
In the corresponding problem, due to longitudinal inertia, the 
longitudinal buckling parameter for a semichord of a uniform 
wing due to a uniform mass distribution and a concentrated mass 
at the leading edge is* 4 

k, = (dc?/4EI)(Mw + 6.33M7) = 15.68 (4) 


Because most wings have low mass/stiffness ratios, this problem 
should not be critical for practical wing designs. 

Although longitudinal inertia effects would appear to be 
negligible for practical wing designs, these effects might be 
more significant for slender fuselages. Only a simple solution 
is presented here, using the assumptions that only one body re- 
gion is flexible, all rigid incidences and deformations are small, 
drag of body acts along the thrust axis, and changes in aero- 
dynamic coefficients, etc., due to deformation are not explicitly 
considered. The notation used is explained by Fig. 1. 

The equations of motion for the configuration simplify to 


Bq = Lifer — gela — Lila + Oe; + ge) — (T — D)Z (5) 
Ma = Ina + Li(a + 6) + (T — D)o (6) 
M.d = (T — D) (7) 


where @ is the angular deformation between the two rigid sections 
land 2. 
The deformation equation is 


Ke6 = Li(a + Oe, + 


4 
f [mig(x + ge) — ma + md(a + @)|xdx (8 
0 


Static instability would occur at some positive value of 6/a and 
Kg at which g = 0. Therefore, by rearranging Eqs. (5)-(8 
and using the result 

Z = Mig.id/M. (9 
the critical joint stiffness Kg is obtained. 


RoL.s. = (L,L2/M,)[Moesee1 - M,s,e2] + 
(Migid/M.)[M-Lose + go MLe — MeL,)| (10 


To assess the effect of longitudinal acceleration d, consider the 
simple example in which Z2/L,; = M2/M, = land e = gi, ie, 
s: = 0. Therefore, the second expression on the right hand side 


of Eq. (10) equals the first when 
d = L,/2M, = L,./2M. (11 


Thus, it is seen that longitudinal acceleration has the same effect, 
in this example, as the lateral aerodynamic and inertia forces 
when the longitudinal acceleration is one half the possible lateral 
acceleration of the rigid body for one radian incidence. 
The condition from Eq. (10) that longitudinal inertia effects 
be negligible is given by 
M,Les2 = go( Mol, — M,Le) (12) 


or 
M,Lre2 = MeL cge (13 
Also, the further condition from Eq. (10) which, with Eq. (13), 
makes the critical joint stiffness zero is 
Mose, = Msye2 (14 


but the satisfaction of both Eqs. (13) and (14) may be difficult 
in practice. 
CONCLUSIONS 

Longitudinal inertia forces should have no significant effect 
on static aeroelastic problems of wings but can have a pro- 
nounced effect on similar problems for aircraft bodies. This 
is because the latter have larger values than the former of the 
ratios mass/stiffness and mass/lift. 

Lateral inertia forces can modify the results of static aero- 
elastic analyses for unconventional wings and should then be 
considered, as they are already in the problems of entire bodies 
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Formulas for the Determination of Thermal 
Stresses in Rings 


Marvin Forray 
Design Specialist, Structures, Applied Research and Development, 
Republic Aviation Corporation, Farmingdale, L.I., N.Y. 


October 21, 1959 


ete PLANE-STRESS, linear elastic solution of a traction-free 
circular ring of inner radius a and outer radius b subjected to 
a general two-dimensional temperature distribution was con- 
sidered in a previous note.! Design equations will now be pre- 
sented for the determination of the polar-coordinate stress com- 
ponents corresponding to 
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READERS’ 


EaT = Ax, nr*® cos n6, 
(K = 0, 1, 2, 


Case A: K+2-—n #0. The stress components are given by 


2261.2...) @ 


Orr = (ao/r?2) + 2ho + [2bir — (2a,’/r*)] cos 6 + 
© 
. | 
> [((m — n?)ayr"™~? + (n + 2 — n?)bar® — 
n=2 | 
(n + n?)an’r®~2 + (2 — n — n*)b,'r-"] cos n6 — 
{[Ax.n(K + 2 — n®)r* cos n6] + 


(K +2+n\(K+2—-n)}} 


org = [2hir — 2(a1’/r*)] sin 8 + > [+n(n — 1)ayr"~? + | 
n=2,3 . 
n(n + l)bar™® — n(n + 1)a,’r-"-2 + n(1 — n) X 


ba’r—"] sin n@ — | [n(K + 1)Ax nr*® sin n6] + 


(K+2+n\(K+2-—n)]} 


oe = (—ao/r?2) + 2b0 + [6bir + (2a;’/r*)] cos 86 + 
[n(n — 1)adnr”™~2 + (nm + 2)(n + L)bar™ + 
em2.3... 
6 ! 
n(n + 1) an’r~"-2 + (2 — n)(1 — n)bn’r~"] cos nO — 


Ax. n»{((K + 2)(K + 1)r* cos n6]/[((K +2+n) xX 
(K +2 —-n)j} 
(2) 


where the coefficients do, bo, b;, a1’, Gn, bn, Gn’, and b,,’ are deter 


mined from the boundary conditions 


Orr = ong = Oonr = a,7 = 5 (3) 
Case Ay: n = 0(K + 2 — n # 0 is automatically satisfied. ) 


This is the axisymmetric case. Substitution of Eqs. (2) into 


(3) yields 


a [—Ax,0o/(K + 2)]} [(6* — a*)a*b?] /(b?—a?)} 
2bo = [Ax,0o/(K + 2)] [(b* *2 — aX *2)/(b? — a?)] (4) 
e' = > = 0, a, = 6, = a,’ = 5,’ = Ofore = 2.3, 4.... 
If T = T(r/b)*, the stress components in nondimensional 


form become 
Or/EaTo = —[(a/b)*/(K + 2)] X 
((b2/r?){ [(b/a)* — 1]/((b/a)? — 1} - 
t ((b/a)**2 — 1]/[(b/a)? — 1]} + (r/a)¥) 
079/EaT, = 0 (5) 
o9/EaTo = —[(a/b)*/(K + 2)] X 
(( —b?/r?){ [(b/a)* — 1]/[(b/a)? — 1]} — | 


} [(b/a)**+2 — 1]/[(b/a)? — 1]} + (K+ 17 a)*) 


Case Ao? n= 1. (K + 2 — 2 # 0 is automatically satisfied.) 
If T = T)(r/b)* cos 6, the stress components are 


6rr/EaTo = {1 (K + 3)] ; (1 — (a/b)*® +3] {1 —(a b)s] fre 


| [(a/b)* +3 — (a/b)4]/[1 — (a/b)*](r/b)3} —(r b)*) cos 6 


o76/EaTy = (1/(K +3)](f{ [1 — (a/b)* *3]/[1 — (a/b)4] }7% + | 
i (a/b) +8 — (a/b)*]/[1 — (a/b)*](r/b)3} — (r/b)*) sin 0 
o9/EaTy = [1/(K + 3)] G3[1 — (a/b)**]/f1 — 
(a/b)*}}"” + { [(a/b)* — (a/b)* +3)/[1 — (a/d)*](r/b)3} — 
(K + 2\(r b)*) cos 6 
(6) 


Case Az: n>1. (K+2-—n #0.) With the substitution 
y = b"-Kb,/Ax, a 


zs = 5-*-*k 2an'/A K,n w= bo” ~Kp,,’ ‘Ag. a 


x = 6*-*-%q,/Axk, 


there results 


FORUM 239 


Orr/EaTy = {x(n — n*)(r/b)"-? + y(n + 2 — n®) X 
(r/b)" + 2(—n — n*)(r/b) "-? + w(2 — n — n*) X 


(r/b)-" — [((K + 2 — n*)(r/b)¥]/((K + 24+) X 


(K + 2 — n)]} cos né 
ore/EaTo = {x(n)\(n — 1)(r/b)"-2? + y(n)(n + 1)(r/b)* — 
a(n)(n + 1)(r/b)-"-2 + w(n)(1 — n)(r/b)-" — 


[n(K + 1)(r/b)*]/[((K + 2+ n)\(K + 2 — n)]}} sin n@ 


o99/EaTy = {x(n)(n — 1)(r/b)"-? + y(n + 2)(n + 1) X 
(r/b)" + 2(n)(n + 1)(r/b)-*-? + w(2 — nl — 1m) X 
(r/b)-" — [((K + 2)(K + 1)(r/b)*)]/[((K + 2 + 1) X 

(K + 2 — n)]} cos n6 
(7) 


where x, y, 2, and w are determined from 
(n — n?)(a/b)"-*-2x + (nm + 2 — n®)(a/b)n Ky 4+ 
(—n — n?)(a/b)-"™-*-2g + (2 — nm — n*)(a/b)” w= 


{(K +2 — n*)/[((K +2+n\(K +2 —n)}} 


(n — n?)x + (n + 2 — n*)y + (—n — n?*)z + 
(2—n—n*)w = {(K+2—n*)/(K+2+n)xX | 
(K +2-—n)}} 
n(n — 1)(a/b)"-*-2x + n(n + 1)(a/b)"-*y — 
n(n + 1)(a/b)-"®-*-22 + n(1 — n)(a/b)-"-*w | 
{n(K + 1)/[(K +2+ n\(K + 2 —-—n)}} 
n(n — 1)x + n(n + 1)y — n(n + 1)2 + nl — n)w 


{n(K + 1)/[(K +24+ n\(K +2 —n)]} 


(8) 
Case B: K+2-—n=0. The components of stress are 
Orr/EaTo = (x"( n — n*)(r/b)"~2? + y(n + 2 — nn?) X 
(r/b)" + 2'(—n — n*)(b/r)"*? + w(2 — n — n*) X 
(r/b)—™ — (r/b)"-2 { [C1 — n)/2] log (r/a) + 


(1/4)[(1/n) + 1}}) cos no 


orep/EalTy = (x'(n? — n)(r/b)"-2 + y’(n)(n + 1)(r/b)" + 
2'(—n — n*)(b/r)"*? + w'n(1 — n)(r/b)-" — (r/b)"-? X 


{ [(n — 1)/2] loge(r/a) + (1/4)[(1/n) + 1)}) sin Ng 


ogo/EaTy = (x'(n)(n — 1)(r/b)"-2 + y(n + 2)(n + 1) X 
(r/b)" + 2'n(n + 1)(b/r)"*? + w'(n — 2)(n —1)X | 
(r/b)~-"™ — (r/b)" 2fI(n — 1)/2] log(r/a) + 


(1 4) [3 — (l ‘n)}}) cos nO 


(9) 
where x’, y’, z’, and w’ are determined from 
(n — n*)x’ + (nm + 2 — n*)(a/b)*y’ — (n + n?) X 
(a/b)—2"2’ + (2 — n — n?)(a/b)?-*w’ = 
(1/4) [(1/m) + 1] 
(n — n*)x’ + (n + 2 — n?)y’ — (nm + n?)2’ + 
(2 —n — n*)w’ = }[(1 — n) log&b/a)] /2} + 
(1/4) [(1/n) + 1] 
(10) 
(n? — n)x’ + n(n + 1)(a/b)*y’ — 
(n + n?)(a/b)~"2z’ + n(1 — n)(a/b)?-*"w’ = 
(1/4) [(i/n) + 1] 
(n? — n)x’ + n(n + 1)y’ — (m + m?)s’ + 
n(1 — n)w’ = { [(m — 1) logb/a)]/2} + 


(1/4)[(1/n) + 1 
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REMARKS 


(1) All stress components are identically zero when K = 2, 
i.e., ZT = To(r/b)" cos nO because r” cos n6 is the real part of the 
analytic function 2" and must satisfy Laplace’s equation V?7' = 0. 

(2) When the temperature distribution is expressed as a series 
of terms of the form (1), the stress components corresponding to 
the individual terms may be superposed to obtaia the resultant. 

(3) The prescribed temperature can always be expressed as the 
sum of symmetrical and antisymmetrical terms in 6. The ex- 
pression (1) represents the symmetrical component, while the 
antisymmetrical components will be made up of terms of the 
form EaT = > Brurr® sin né. 

Kn 

The above formulas remain valid if cos 2@ and sin 76, respec- 

tively, are replaced by sin v6 and —cos n8@, respectively. 
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Stress Functions for Plates Bounded 
by Piecewise Analytic Curves 
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N A RECENT NOTE, Balmer! has shown that the stress function 
for a plate of polygonal shape fulfilling the stress-free bound- 
ary conditions can be derived from the generalization of the 
standard expression for the Airy stress function for rectangular 
plates. It appears that this generalization can be carried out 
further to include planforms which are bounded by piecewise 
analytic curves, Fig. 1. 
Suppose the plate is bounded by m analytic curves which can 
be expressed as 
g(x,y) = 0, ¢=1,2,....™m (1) 
The stress function whichful fills the stress-free boundary condi- 
tions on each curve can then be expressed in the following form 


m 


o(x, y) = II [gi(x, y)]2P(x, y) (2) 
t=1 


where P(x, y) is an arbitrary function which is nonvanishing over 
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Fic. 1. Arbitrary planform bounded by analytic curves. 
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Fic. 2. Planform bounded by an ellipse and a straight line. 


planform. The stress-free boundary conditions require that 


0° 


0° 
On ;dl; 


( Nne Mi 


where (JN,,); and (V,,); are the normal and shear stress resultants 
and n; and ¢; are the normal and tangent on the ith boundary, 
respectively. 

In order to show that Eq. (2) satisfies the condition of Eq. (3), 
it suffices to consider only one edge. For the ith edge, we write 


o(x, y) = [gi(x, y)]2Q(x, ¥) (4) 


where Q(x, y) includes P(x, y) and the remaining functions. The 
second derivatives in directions of m and tare 
Nn = o,te = $,2r2(OX/Ot)? + ,yy(Oy/OL)? + 2,2y(Ox/Ot)(Oy/OdE) 
Ni = Onn = %,22(Ox/On)? + 
,yy(Oy/On)? + 2¢,7,(Ox/On)(Oy/On) 
—Nnt = Ont = ,22(Ox/Ot)(Ox/On) + o,yy(OV/Ot)(Oy/On) + 
$,2y(Ox/Ot) (Oy/In) + (Oy/Ot) (Ox/On) (5) 


Differentiating ¢ of Eq. (4) with respect to x and y and noting 


that gi(x, vy) = Oon curve 7, one obtains 


d,rxr = 2Q[0gi(x, y)/Ox]? 
P.yy = 2Q[dgi(x, y)/Oy]? 
dry = 2Q} [Ogi(x, y)/Ox][dgi(x, v)/dy] } (6) 


Substituting Eq. (6) into Eq. (5), one obtains 


(Ni)i = 2Qf [Ogi(x, y)/Ox] (Ox/dt:) + [Ogi(x, ¥)/Oy](Oy/ dt:)} 
= 2Q[0g;/Ot;]? 

(Ni)i = 2Q[0g;/On;]? 

(Nnt)i = —2Q[(0g:/0n;)(O0g:/Ot;)] 

Since gi(x, y) = 0 on curve i, it is obvious that 0g;/dt; vanishes. 

Consequently the free-stress boundary conditions 


(Nn) = 0, (Nat)i = 0 


are satisfied. 

As a simple example, the stress function for a plate shown in 
Fig. 2 can be expressed as: 

o(x, y) = [y + (c/d)x + c]*(b2x?2 + a?y? — a2b?)?P(x, y) (8) 

Aside from the stress-free boundaries, there can be boundaries 
on which either the stresses or the displacements are prescribed. 
These boundary conditions can be satisfied on the average 
through the use of a variational principle. However, for the 
case of prescribed stresses, it is sometimes possible to satisfy 
these conditions by proper choice of the function P(x, y). 
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THE B, F. GOODRICH COMPANY 

GOODYEAR AIRCRAFT CORPORATION 

GRUMMAN AIRCRAFT ENGINEERING CORPORATION 


HYDRO-AIRE CO. 
INSURANCE COMPANY OF NORTH AMERICA COM- 
PANIES 


INTERNATIONAL BUSINESS MACHINES CORPORATION 
THE INTERNATIONAL NICKEL COMPANY, INC. 


ITT ——— DIVISION OF INTERNATIONAL 
TELEPHONE AND TELEGRAPH CORPORATION 


JANITROL AIRCRAFT, A DIVISION OF MIDLAND- 
ROSS CORPORATION 


JOHNS-MANVILLE SALES CORPORATION 

WALTER KIDDE & COMPANY, INC. 

KOLLSMAN INSTRUMENT CORPORATION 

LAVELLE AIRCRAFT CORPORATION 

LEAR, INCORPORATED 

C. W. LEMMERMAN, INC, 

LIBRASCOPE DIVISION OF GENERAL PRECISION, INC. 
THE LIQUIDOMETER CORPORATION 

LOCKHEED AIRCRAFT CORPORATION 


LOEWY-HYDROPRESS DIVISION OF BALDWIN-LIMA- 
HAMILTON CORPORATION 


THE MARQUARDT CORPORATION 

THE MARTIN COMPANY 

McDONNELL AIRCRAFT CORPORATION 

MELETRON CORPORATION 
MINNEAPOLIS-HONEYWELL REGULATOR COMPANY 
NATIONAL CREDIT OFFICE, INC. 

NORTH AMERICAN AVIATION, INC. 

NORTHROP CORPORATION 


NUCLEAR DEVELOPMENT CORPORATION OF AMERICA” 
PAN AMERICAN WORLD AIRWAYS, INC. 
THE RALPH M. PARSONS COMPANY 

pas *~ pepe DIVISION, BORG-WARNER CORPOReS 


PHILLIPS PETROLEUM COMPANY 

PIASECKI AIRCRAFT CORPORATION 

RADIO CORPORATION OF AMERICA 
ASTRO-ELECTRONIC PRODUCTS DIVISION 

DEFENSE ELECTRONIC PRODUCTS 


RAMO-WOOLDRIDGE DIVISION, THOMPSON 
WOOLDRIDGE INC. 


REPUBLIC AVIATION CORPORATION 
ROHR AIRCRAFT CORPORATION 
PAUL ROSENBERG ASSOCIATES 
RYAN AERONAUTICAL COMPANY 
SANDBERG-SERRELL CORPORATION 


SHAFER BEARING DIVISION, CHAIN BELT COMPAR ; 


SHELL OIL COMPANY 

SIMMONDS AEROCESSORIES, INC. 

SOCONY MOBIL OIL COMPANY, INC. 
SOLAR AIRCRAFT COMPANY 

SOUTHWEST PRODUCTS CO. 

SPACE TECHNOLOGY LABORATORIES, INC. 


R. DIXON SPEAS ASSOCIATES 


SPERRY GYROSCOPE COMPANY DIVISION OF SPERRy 
RAND CORPORATION 


STANDARD OIL COMPANY OF CALIFORNIA 
STANDARD OIL COMPANY (INDIANA) 


so es ae CORPORATION 
ORD MANUFACTURING DIVISION 


STANLEY AVIATION CORPORATION 
THIOKOL CHEMICAL CORPORATION 
THOMPSON RAMO WOOLDRIDGE INC, 
TRANS WORLD AIRLINES, INC. 

TURBO PRODUCTS, INC. 

UNION CARBIDE CORPORATION 

UNITED AIR LINES, INC. 

UNITED AIRCRAFT CORPORATION 
UNITED STATES AVIATION UNDERWRITERS, INC. 
VERTOL AIRCRAFT CORPORATION 
VITRO CORPORATION OF AMERICA 
WESTERN GEAR CORPORATION 
WESTINGHOUSE ELECTRIC CORPORATION 
WYMAN-GORDON COMPANY 

YOUNG RADIATOR COMPANY 
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